Evaluating Results from the Relativistic Heavy Ion Collider 
with Perturbative QCD and Hydrodynamics 



R. J. Fries, 1 '2 C. Nonaka^ 

^Cyclotron Institute, Texas A&M University, College Station TX, USA 
^RIKEN/BNL Research Center, Brookhaven National Laboratory, Upton NY, USA 
"^Department of Physics, Nagoya University, Nagoya, Japan 

December 15, 2010 



Abstract 



We review the basic concepts of perturbative quantum chromo dynamics (QCD) and relativistic 
hydrodynamics, and their apphcations to hadron production in high energy nuclear collisions. We 
discuss results from the Relativistic Heavy Ion Collider (RHIC) in light of these theoretical ap- 
proaches. Perturbative QCD and hydrodynamics together explain a large amount of experimental 
data gathered during the first decade of RHIC running, although some questions remain open. 
We focus primarily on practical aspects of the calculations, covering basic topics like perturbation 
theory, initial state nuclear effects, jet quenching models, ideal hydrodynamics, dissipative cor- 
rections, freeze-out and initial conditions. We conclude by comparing key results from RHIC to 
calculations. 
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1 Introduction 



The Relativistic Heavy Ion Collider (RHIC) at Brookhaven National Laboratory started operations 
about a decade ago. The amount of data collected and the quality of the data have been outstanding. 
Besides a successful proton-proton and proton-nucleus program, RHIC has mostly provided data on 
nuclear collisions, from a few GeV center of mass energy up to 200 GeV per nucleon-nucleon pair. We 
have strong evidence that the central goal of the RHIC program, the discovery of quark gluon plasma 
(QGP), a deconfined state of nuclear matter, has been achieved. In order to draw this conclusion a 
wide variety of observables have been weighed against theoretical expectations and we will discuss a 
few of those in this article. Some key experimental discoveries at RHIC over the past decade were 
(i) the extremely strong jet quenching, many times that of ordinary nuclear matter [1] (ii) the 
very large elliptic flow of the fireball that confirms collective behavior at energy densities larger than 
expected at the phase transition [3]; (iii) the surprising quark number scaling of elliptic flow that seems 
to indicate that the collective flow is carried by quarks [H [5] (see [6] for an attempt of an alternative 
explanation); and (iv) direct photon measurements that suggest large initial temperatures [7]. Even 
before the partonic nature of the fireball could be established there was mounting evidence that the hot 
matter at RHIC was not behaving like a weakly interacting gas, but rather like a strongly interacting 
liquid. This has led to the conjecture that quark gluon plasma is a nearly perfect liquid [8], at least 
close to the phase transition temperature. Conservative estimates for the initial energy density in the 
center of head-on collisions at top RHIC energy find a lower bound ~ 3 GeV/fm^, which is above the 
estimated critical energy density [3]. 

Perturbative quantum chromodynamics (pQCD) and relativistic hydrodynamics have been two im- 
portant tools to understand and interpret RHIC data. It was found that the bulk of the produced 
particles at RHIC (for transverse momenta Pt smaller than ^ 2 GeV/c) show signatures of collective 
behavior. The mean-free path of particles seems to be small enough for the dynamics to be described 
by relativistic fluid dynamics. This was a non-trivial finding since hydrodynamic descriptions for lower 
energy nuclear collisions routinely overestimated the amount of collectivity. While hydrodynamic mod- 
eling 10 years ago was still rough, based on (2+l)-dimensional ideal fluid dynamics with simple initial 
conditions and freeze-out, there has been an amazing amount of progress since then by going to full 
(3+l)-dimensional modeling, taking into account dissipative corrections, fine-tuning of initial condi- 
tions all the way to event-by-event calculations, and a deeper understanding of the hadronic phase with 
separate chemical and thermal freeze-outs, and through the advent of hybrid hydro-|-cascade models. 
We will highlight many of these improvements in this article. The progress has enabled hydrodynamic 
models — and the entire RHIC program — to enter a phase in which quantitative measurements are 
finally close. Prime candidates for such quantitative measurements are the equation of state of hot 
QCD, including the order of the phase transition between hadronic matter and QGP and the existence 
and location of a critical point, and the shear viscosities of these phases. The measurement of other 
bulk transport coefficients, like the bulk viscosity and relaxation times, are in principle possible but 
remain elusive for now. We will discuss the status and potential problems of such measurements. 

Hydrodynamics describes the bulk of the particles in a collision (more than 98% of them). The tail 
of the particle P-r-spectra in nuclear collisions, which clearly contain particles that have not thermalized, 
should not simply be disregarded. In fact it was proposed a long time ago that they can serve as "hard 
probes" of the bulk matter created. In elementary p + p oi p + p collisions hadrons with transverse 
momenta of 5 GeV/c or more are created through a single hard scattering of two partons within the 
wave functions of the colliding hadrons, which then fragment in the vacuum away from the collision into 
coUimated bunches of hadrons, called jets. This entire process can be calculated in perturbative QCD 
due to the large momentum transfer involved, while the unavoidable non-perturbative contributions can 
be treated in a controlled way through a formalism called collinear factorization. Perturbative QCD 
based on collinear factorization has been a great success story in elementary collisions |9j. Hard initial 
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scatterings of partons from the initial nuclear wave functions should proceed in a way very similar to 
elementary collisions, with the understanding that the wave functions of free nucleons and those in nuclei 
might differ somewhat. However, the big difference arises in the final state, when an outgoing parton 
or jet finds itself embedded in a fireball of hot and dense quark gluon plasma. Clearly we expect those 
partons to rescatter and lose energy through elastic collisions or bremsstrahlung. The final state effects 
on high-Pr hadrons and jets should encode valuable information about the QGP phase. The most 
prominent example is the transport coefficient q that parameterizes the average momentum transfer 
per unit path length to a fast parton in the medium. The so-called LPM effect, coming from the finite 
formation time of induced radiation, leads to a signature quadratic dependence on the thickness of the 
medium. We will focus our attention here on the leading particle description which has received the 
most attention by theoreticians and is the most relevant effect for observables measured so far. Despite 
this restriction to the apparently simplest problem we will see that we do not yet have a consistent 
description of this problem. We will not deal in detail with more comprehensive approaches that follow 
full jet showers in the medium. 

In this review we want to lay out the basic concepts of both perturbative QCD and relativistic 
hydrodynamics and their applications to hadron observables in nuclear collisions at RHIC. This will 
enable us to discuss some important results from RHIC and to draw conclusions. The article is organized 
as follows. In Section [2] we review the fundamentals of collinear factorization, parton distributions and 
fragmentation functions and simple pQCD cross section computations. We will then see how these 
processes change in a nuclear environment leading us to nuclear shadowing and the Cronin effect. We 
then proceed to discuss final state energy loss and the LPM effect. We focus on four common models 
of leading parton energy loss. Finally we give a quick overview of photon production in heavy ion 
collisions. In Section [3] we present the basic concepts of both ideal and viscous hydrodynamics, and 
quickly comment on possible numerical implementations. Then we connect hydrodynamics to the bigger 
picture of heavy ion collisions and discuss necessary details like the equation of state, initial conditions 
and freeze-out procedures. We also briefly touch upon quark recombination. In Section H] we discuss 
data from RHIC in light of the previous two sections. We present single particle spectra for light 
hadrons, azimuthal asymmetries, hadron correlations, and photons and their correlations. We have 
omitted heavy quarks and dileptons in this review which are very interesting topics in their own right 
but would have significantly increased the size of this article. Section |5] contains our conclusions and 
summary. Along the way we try to emphasize practical applications of theory over technical derivations. 
We hope that this article serves as a useful guide for the practitioner. 

2 Particle Production in Perturbative Quantum Chromody- 
namics 

Perturbation theory is a well established tool to deal with interacting quantum field theories. In 
quantum electrodynamics (QED) it has produced some of the most accurate predictions confirmed by 
experimental data. The basic concept is an expansion of observables in powers of the coupling constant 
g of the theory if (7 ^ 1. Naturally, this method becomes unreliable if g is too large. Unfortunately, in 
quantum chromo dynamics the strong coupling as = g^ j (47r) grows logarithmically as the momentum 
transfer squared decreases. This behavior immediately raises serious questions about the usefulness 
of perturbation theory in any realistic situation. Weak coupling methods should work in the asymptotic 
limit — ?► 00. But QCD bound states, hadronization, and the transport properties around the QCD 
phase transition temperature are completely outside the perturbative region. Nevertheless perturbation 
theory in QCD, truncated after the few lowest orders, together with a rigorous factorization program, 
to separate off infrared divergences representing the long distance behavior of QCD, has been shown to 
work down to ~ 1 GeV^ in some applications. Not all processes feature a rigorous and unambiguous 
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factorization, but we have hope that hadron and photon production at large transverse momentum 
Pt in nuclear collisions can be described by perturbative methods. The goal of this program is it to 
use high-Py hadrons and jets as hard probes, whose final state interactions with the bulk of the event 
reveals important information about quark gluon plasma. 

In the first subsection we review some of the basic principles of computing the cross sections or yield 
of high-Py hadrons and jets at hadron colliders. In the second part we discuss modifications expected 
in collisions involving nuclei, focusing in particular on initial state, or cold nuclear matter effects. In the 
third part we address final state effects and parton energy loss for hadrons and jets, which we critique 
and compare further in the fourth part. In the last subsection we briefly discuss the production of 
photons. 

2.1 Factorization in pQCD 

Even though the creation of hadrons at large transverse momentum Pt ^ 1 GeV involves a large 
momentum transfer Q P^, one has to deal with the fact that the initially colliding hadrons Hi and H2, 
and the final hadrons H^, H4, ... are multi-parton states in QCD bound by non-perturbative dynamics. 
Fortunately, for several key processes it has been possible to prove factorization theorems [101 [HI [121 
[T3l [HI [15] , see also [9] for a more didactic introduction. They allow us to separate perturbative and 
non-perturbative processes in a well-defined, systematic way by factorizing all infrared and long-range 
dynamics into universal, well-defined and observable matrix elements. They establish an expansion (in 
powers of l/Q) of the underlying "hard" partonic process. The leading process in l/Q, often called 
leading twist, is usually the one with the fewest possible partons connecting to the long-range part. E.g. 
for single (or di-hadron) production from two hadrons, H1 + H2 — > H^ + X (or H1 + H2 — )■ H3 + H4 + X), 
the leading underlying parton process is that of 2 — )■ 2 scattering of partons a + 6 — ?■ c + d with parton 
a, b, c, (d) being associated with bound states Hi, H2, H^, (H^), resp., see Fig. [H The blobs in Fig. 
[T] represent the association of one parton with its parent hadron. In the initial state these are called 
parton distributions in the final state they are fragmentation functions -Dc/h- 
The factorization of the cross section can schematically be written as 



where aa+b^c+x is the hard partonic cross sectioiu involving a large momentum transfer. Processes 
involving other "associations", most notably those with more partons taken from one hadron are of 
higher twist and suppressed by powers of l/Q ~ 1/Pr- The convolution signs mean that the parton 
momenta connecting blobs and hard cross sections have to be integrated, if not fixed by kinematics. 
This will become more clear when concrete examples are discussed further below. 
A few additional remarks are in order. 

• Here we will only deal with collinear factorization. This is sufficient for processes with a single 
hard scale Q, or even for processes with different scales Qi, Q2 (then resummations are needed) as 
long as all scales are large. So called /cT-factorization is needed for processes with one hard and one 
soft scale and will not be discussed here [TB1[T7|. Practically this means partons can kinematically 
always be treated as collinear with their parent hadrons, which simplifies the momentum integrals 
in the factorization formulas tremendously. 

• At first we will only discuss leading twist processes. With scales of the order of a few GeV this 
is sufficient for single hadrons. However, in the case of nuclei we will see that some higher twist 
processes are enhanced and become important. 

^More precisely this is the cross section modulo some collinear and infrared divergences which have been factorized 
into the parton distributions and fragmentation functions 
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Figure 1: Schematic sketch for the amphtude of the leading factorized processes for pro- 
duction of a single hadron in the collision of two hadrons Hi and H2- The thin lines 
represent single partons connecting the soft (i.e. long distance) processes within the hadrons 
with the hard (i.e. short-distance) scattering. Soft and hard processes are depicted by round 
blobs and boxes respectively. The probability distributions referring to the soft and hard 
processes are defined as the parton distributions / or fragmentation functions and the 
hard partonic cross section aa+b~>c resp. 



• We will also refrain from discussing particle production in the limit of very large center of mass 
energy, when the gluon distribution of hadrons saturates. For large nuclei this limit might be 
reached at RHIC energies and particle production from this Color Glass state could be dominant 
for particles at lower Pt (from scattering of partons with low Bjorken-^ in the initial wave func- 
tion). With the saturation scale Qs for gold nuclei at RHIC energies estimated to be smaller than 
2 GeV [l8j, the high-Py domain should be safely in the region of collinear factorization. We will 
revisit this topic in the section about initial conditions for hydrodynamics. For the most recent 
reviews of the Color Glass Condensate see e.g. [T9l [20] . 

• Collinear factorization has been rigorously proven in very few cases, and even for simple processes 
there are examples where factorization breaks at a certain order in l/Q" [21j. This is particu- 
larly worrisome for collisions of nuclei where multiple scattering and higher twist corrections are 
enhanced. We will always assume factorization of initial states and hard processes here. On the 
other hand the study of final state effects in nuclear collisions is by definition an investigation of 
how factorization and universality are broken for long-distance final states. 



2.1.1 Cross Sections of Partons 

The factorization theorems mentioned above make sure that the underlying hard parton cross sections 
are infrared-safe. They can be calculated in a perturbative expansion. Singularities from radiative 
corrections can be factored off into the long-distance part that is described by parton distributions and 
fragmentation functions. E.g., for our example of single hadron production at large momentum the 
underlying parton cross section can be written as 
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Leading order (parton model) cross sections are easily calculated. For further processing parton 
cross sections a + 6 — )■ c + c? are most easily parameterized in terms of the Lorentz- invariant Mandelstam 
variables s = (Pa + Pb)'^, t = {Pa—Pc)^, u = {pa — PdY where Pa, Pb, etc. are the four-momenta of partons 
a, b, etc. For example for the scattering of two different quark (or antiquark) species q and q' 



dt s2 2m t2 



(3) 



where A^c = 3 is the number of colors, and by definition we have averaged over ingoing spins and colors 
and summed over outgoing spins and colors (we have kept the coupling constant as part of the cross 
section unlike indicated in ([2])). A comprehensive table for production of light partons and photons can 
be found in the review article by Owens [22]. 

Next-to-leading (NLO) calculations of parton production is much more challenging. The basic matrix 
elements can be found in the work by Ellis and Sexton |23]. The one- and two-jet cross sections were 
e.g. worked out in Several numerical codes are available for jet, hadron or photon production 

at NLO accuracy, performing the required phase space integrals and cancellation of singularities. An 
excellent starting point for the interested reader is the PHOX collection by Aurenche and collaborators 
[26]. 

We have to discuss an important point here. From the NLO-level on cross sections with parton 
final states are no longer well defined, i.e. infrared-safe. In fact we can only define cross sections either 
into hadrons or jets, i.e. sprays of hadrons defined by energy in a restricted region of phase space. At 
leading order one can make the convenient identification jet = parton. At NLO two partons can be so 
close together in phase space that they have to be replaced with one jet. 

In nuclear collisions with its emphasis on final state effects the convenient identification becomes a 
necessary simplification that allows for the treatment of energy loss and other effects on the basis of 
single partons. This is also one reason why a large fraction of literature on heavy ion collisions uses 
leading order calculations. Recently, more and more NLO-based calculations have been presented. In 
that case caution is in order if they are combined with final state effects based on a single parton picture. 

The basis for the use of LO cross sections is the fact that for single and double hadron and photon 
Py-spectra LO accuracy yields reasonably good results. It turns out that for collisions of single hadrons 

with a K-factoT that is close to one and only weakly dependent on the momentum P^ of the produced 
hadron [2?]. For convenience K is hence often approximated by a constant. 



2.1.2 Parton Distributions and Fragmentation Functions 

In the schematic factorization formula ([T]) fa/Hi and fb/H2 are parton distribution functions (PDFs) 
which describe the probabilities that partons a, b can be found in hadrons Hi, H2, resp., with given 
momenta. Note that factorization at leading twist provides a very satisfying probabilistic picture (there 
is no interference between amplitude and complex conjugated amplitude of the parton line connecting 
the hard cross section with bound states). Parton distributions are well-defined and gauge invariant 
matrix elements in QCD. They are also universal, i.e. their definition is independent of the particular 
process in which they occur. 

Suppose hadron H is moving with large momentum P along the positive z axis such that P~^ — )■ 
00. We introduce the light cone components of a four-vector p^ as = (p° ±p^)/\/2. The parton 
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Figure 2: Left panel: diagrammatic representation of a quark parton distribution. The two 
"end" points represent the positions of the quark fields in Eq. (|5]). Right panel: diagrams 
with lowest order radiative corrections for the quark parton distribution in light cone gauge. 
Note that there are generally more diagrams in different gauges due to the presence of a 
gauge link. 



distribution for quarks and gluons in a light cone gauge (A"*" = 0) are defined as 

= / ^e-'^^"^- {H{P)\q{y-)^U{^) \H{P)) (5) 

= {H{P)\F^''{y-)Fa.-'{0) \H{P)) (6) 

Here \H{P)) is the suitably normalized single hadron state, (note that an averaging over hadron spins 
is usually silently implied in the notation {H{P) \ . . . \H{P))), and q and F are the operators for quark 
and gluon fields. ^ with < ^ < 1 is the momentum fraction of the parton in the parent hadron. In 
light cone gauge it is straight forward to interpret these matrix elements as quark and gluon counting 
operators. Note again that these matrix elements can not be evaluated perturbatively for hadrons or 
nuclei. The left panel of Fig. |2] shows a diagrammatic representation of a parton distribution function 
in light cone gauge. 

Radiative corrections introduce a weak, logarithmic scale dependence. The first radiative corrections 
are shown in the right panel of Fig. [2] for a light cone gauge. Resummation of these diagrams lead to 
the DGLAP evolution equations which determine the running of parton distributions fa/H{^,fJ') with 
the scale fi [281 ISll [30| . For a parton a they are 



dfa{i,lj) _ as dy ^ ^^^^ 



(7) 



where the set of Pb^aiv) are called the splitting functions. This notion is easily explained with a look 
at the right panel of Fig. [2] whose diagrams represent the q q splitting function which is 

P,^M = Cf ^ + cj^-6{y - 1) (8) 

where = 4/3 is the color factor. Virtual corrections make a contribution aX y = 1 which introduces 
the (5-function term and regularizes the singularity in the first term via the + description: {f{y)/[l — 
7/])+ — >■ [fiy) — /(1)]/[1 — y] for the integral over any function — y]. More details and the full 

set of splitting functions are discussed in [9J. 

While the /^-dependence is hence perturbatively calculable, the ^-dependence can only be extracted 
from fits to data. This relies heavily on data from the "clean" deep-inelastic scattering (DIS) process. 
Very accurate parameterizations including estimates of uncertainties are available for protons and, via 
isospin symmetry, for neutrons in a wide range of about 10~^ < x < 0.5. The most used parameter- 
izations are from the CTEQ [311 [32] and MRST collaborations [331 El]- The Durham data base has 
comprehensive information about PDFs [35] . Parton distributions of nuclei are discussed further below. 
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Fragmentation functions D(./h{z,^) give the reverse probability that hadron H hadronizes from 
parton c in the vacuum with a certain momentum fraction z of the parent parton [36j. Unhke the 
case of parton distributions a complete sum over states can not be removed and hence fragmentation 
functions can not be written as a single forward matrix element. Instead we have 

D,/h{z,^,) = z I ^e-^^^y-/^H{P)\q{y')\0)^UO\q{0)\H{P)) (9) 

for a quark q with large light cone momentum P+. As for parton distributions, vacuum fragmentation 
functions have been parameterized, mostly from hadron production data in + e~ collisions, but 
uncertainties in the fits are appreciable, even for quite common hadrons like protons and kaons. In 
addition, some sets are not isospin-separated, i.e. they only parameterize processes like u + u — )■ vr^ + vr". 
This uncertainty in the theoretical baseline makes the search for nuclear effects more challenging. The 
most widely used parameterizations in heavy ion physics are the sets by Kniehl et al. (KKP) [37] and 
Albino et al. (AKK) [38], Hirai et al. (HKNS) [39] and deFlorian et al. (DSS) [iQlST]. The latter ones 
include iso-spin separation, partially also including data from p + p collisions. 



2.1.3 Factorized Cross Sections 

In this subsection we will summarize some often used factorization formulas for hadron or jet production 
at leading order. They can be used together with the list of leading order parton cross sections in [22] 
and the parton distributions and fragmentation functions referenced above to make estimates for rates 
of hadron and jet production. Our starting point is the differential production cross section of two 
partons c and d from two hadrons A and B 

dOAB^cd = fa/ A (^g) fb/B {jb) d (^ab—^cd ■ 

(10) 

a,b 

We need to introduce some notation for the kinematics. Let the momenta of the parent hadrons be 
Pa and Pb in positive and negative direction, resp., along the z-axis in the center of mass frame of the 
hadrons. We assume that P = P^ = P^ is larger than any relevant masses. Note that the kinematics 
in ( ITOj) is fixed at leading order with = pt / P and C,b = Pb/P, resp. where pa, Pb, Pc and pd are the 
momenta of the four partons. 

One can easily deduce the cross section for a di-jet event with final rapidities yc and yd and transverse 
momentum pt (the transverse momenta of c and d are equal and opposite), 

dcTAB^cd j: j: /^\^daab^cd x 
2_^Uja/A\ka)kbJb/B\Sb) -7, \^^) 



27ipTdpTdycdyd , ' n dt 

a,b 

where the momentum fractions are fixed to be 

U = — (e^^ + e^^) , 6 = — (e-^= + e"^^) . (12) 
s s 

Note that all Mandelstam variables s, t, u are defined on the level of partons a, 6, c, d. In particular 

t = -iaPT^e-y^ (13) 

where 5* = (P4 + PbY — is the total center of mass energy squared of the two hadrons. 

For single jet events we have to integrate one of the final parton momenta. This introduces effectively 
one non-trivial integral which is usually rewritten as an integral over one of the initial parton momentum 
fractions. For a jet with transverse momentum px and rapidity y we find 

daAB^c+X \^ JC f (c\f ^a^b dcTab^cd 

o 1 — T' ^ 1^ I d^aJa/A{U)Jb/B[^b)-7r^ T — - — -J- — (14j 

2'KpTdpTdy Ti2ia-iTey dt 



9 



where is fixed to 



and the integration boundary (to keep ^j, < 1 for fixed px) is 

= 2-iTe-y' ^^^^ 

We have introduced the useful scahng variable = 2pt/VS- 

In order to arrive at cross sections for hadrons we have to multiply the differential cross section 
for partons with the corresponding fragmentation functions. The resulting phase space integrals are 
often shifted to be integrals over the initial parton momentum fractions C,a and C,b- For applications in 
heavy ion physics we rather adopt a different way that keeps factorizability between the fragmentation 
functions on one hand and the parton cross section plus parton distributions on the other hand explicit. 
For two hadrons C and D with momenta Ptc and Ptc "we can write 

dNnn sr^ dzr. dzd 



2tt PTcdPrcPTDdPTDdycdyD 



E 



X 7. , , , Dc/c{zc)D,ij,{z,) . (17) 

2'KpTcdpTcPTddpTddycdyD 

This is a very general formula that connects a distribution function of partons c, d with momenta 
Pc = Pc/zc, Pd = Pn/zd, resp. in the final state to hadrons C, D. Of course the applicability of 
this formula still requires the collinear fragmentation picture to hold in this much generalized setting. 
Nevertheless, derivatives from Eq. (1171) are often used to model final state interactions for hadron 
production in nuclear collisions. There might be kinematic constraints that lead to lower bounds on 
the integrals over Zc and z^, whose exact specification will depend on the distribution N^d of partons. 

For completeness and further clarification let us discuss the more familiar special case of hadron 
production in a regime where final state interactions can be neglected, e.g. m p + p collisions. The 
formula above holds also for cross sections, Ncd — ^ o'ab-^cd, and the partonic cross section is given 
by ( |TT1) times an obvious phase space factor l/prd^iPTc — Pxd)- This factor can be used to cancel the 
integral over Zd to lead to 



daA 



B-^CD 



2tx PTcdPrcPTDdPTDdycdyD 

1 NT^ Z'^""'' dZr ^ „ . . . . „ . . . 1 „ , . „ , .dan 



PtcPtd , , ■ 

a,b,c,d " ^™ 



/ ^^afa/Ai^a)^bfb/Bi^b)-Dc/c{Zc)Dd/DiZd)^^jp^ 



where 



Zd = Zc— — , 2;max = miu <^ 1, — — ^ , z^i^i = ^=-max {coshyc, coshyz)} . (19) 

Note that and are given by f fT2|) with ptc = Ptc/zc and pTd = Pro/zd- 

For single hadron production we can provide a similar general formula for fragmentation from a 
distribution of partons Nc, 

dNc \- f' dz , dNc 

1^ l^D./c{z) (20) 



27iPTdPTdy ^ z^ 27rpTdpTdy 
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It can be applied if collinear fragmentation is the correct description of hadronization of an ensemble 
of partons and obviously p = P/z. The special case of single hadron production in collisions with 
negligible final state interactions gives the formula 



d(TAB^C+X 

2TiPTdPTdy 




TT 2^a - ^TCy dt 



^a^b daab^cd 



(21) 



where 



2Pt 



cosh y 



(22) 



mm 



and the other kinematic variables can be inferred from f lTSj) and f lT6|) . 

For an alternative way of handling the phase space integrals in terms of hadron production see [22] . 
Let us point out once more that the convenient identification of single partons and jets is only valid in 
the context of leading order calculations. 

2.1.4 Photons 

In principle, photons with high transverse momentum Pt can be treated in a fashion very similar to 
hadrons or jets. We usually do not consider photons from decays of hadrons (predominantly vr'') long 
after the collision. After subtracting those decay photons we are left with the "direct" photons produced 
in the collision. Photon yields produced directly in the hard process can be calculated via Eq. f fT^ 
together with the corresponding parton level processes. At leading order, the cross sections of the 
annihilation and Compton diagrams, q + q ^ + g and (? + g — )■ 7 + g, resp. can be found in the work 
by Owens [22]. 

Another way to produce direct photons is as bremsstrahlung in hard process like g + g ^ q + q. 
One of the outgoing quarks can radiate a collinear photon while fragmenting. This process can be 
described by photon fragmentation functions [2^ Eq. f l20p together with the usual set of hard 
parton processes and photon fragmentation functions are used to compute this contribution. At next- 
to-leading order, bremsstrahlung and hard photon radiation in the final state have to be calculated 
in a consistent scheme to separate large angle and collinear photon radiation [13| HI]. NLO direct 
photon calculations have had some difficulties in the past to describe all aspects of photon production 
in hadronic collisions [15]. The theoretical understanding has been improved in recent years by the use 
of various resummation techniques [l6l H?]. The PHOX codes can also be used for NLO calculations of 
direct photon production. 

Photon-hadron and photon-jet pair production is a particularly hot topic in heavy ion collisions as 
we will discuss in more detail later. Their yields with both initial hard and bremsstrahlung photons 
can also be calculated in a straight forward way from the factorization formulas in the last subsection. 
Fragmented photons can in principle be distinguished from prompt hard photons since the latter are 
not accompanied by hadrons close by in phase space while the former are usually part of a jet cone. 
Experimentally, isolation cuts for photons can help to suppress the bremsstrahlung contribution and 
give access to more detailed information. 

In nuclear collisions there are additional sources of direct photons. We will discuss jet conversion 
into photons in the subsection about final state interactions. There is also thermal radiation from the 
hot hadronic matter, and, if energy densities are large enough, from the partonic QGP phase. In fact 
the latter is one of the key observables that we would like to study at RHIC, since the thermal photon 
spectrum can work as a direct (though time- and space-averaged) thermometer of the quark gluon 
plasma. To compute the photon spectrum the time-evolution of the QGP fireball has to be folded with 
rates as a function of the local temperature and chemical potential. The time evolution is naturally 
done through a hydrodynamic model as discussed in Sec. [3] of this review. The rates can be calculated 
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perturbatively if the temperature T is large enough to render the strong couphng constant small, 
g <^ 1. Although there is some doubt whether the maximum temperatures at RHIC (T < 450 MeV) 
are sufficient to warrant a perturbative description the perturbative computation of thermal photon 
rates has been a sustained effort over many years [HI HU |50]. The complete leading order results have 
been given by Arnold, Moore and Yaffe in [5T| 152] . Additional photon radiation could be emitted in the 
pre-equilibrium phase. In particular, the fact that quarks and gluons are not in chemical equilibrium 
early on could affect photon rates and would not be mimicked well by hydro codes initialized at very 
early times. Estimates for pre-equilibrium photon yields in a transport approach can be found in [53] . 

2.2 Nuclear Collisions: Initial State Effects 

Perturbative techniques and factorization have first been developed for scattering involving individual 
hadrons. The basic principles should be valid if one or both of the scattering partners are bound in 
a nucleus. In fact, the small binding energies and slow relative motion of nucleons should not have 
large impact on scattering at large momentum transfer. However, the larger volumes filled with nuclear 
matter surrounding the point-like hard interaction should potentially lead to rescattering both in the 
initial and final state. In this subsection we will discuss the most relevant nuclear effects in the initial 
state. 

2.2.1 Shadowing and Nuclear Parton Distributions 

Individual nucleons are clearly distinguishable building blocks of nuclei. Hence we expect parton distri- 
butions in a nucleus with Z protons and A — Z neutrons to be well represented by a linear superposition 
of parton distributions of the individual protons p and neutrons n, 

The remaining non-trivial nuclear modification Ra was expected to be small until it was found by the 
EMC collaboration that deep-inelastic scattering off nuclei leads to sizable differences between free and 
bound nucleons [M]- Note that nuclear parton distributions are usually normalized to one nucleon. 

Despite considerably larger uncertainties compared to free nucleon parton distributions we can now 
identify four distinct regions of behavior in the momentum fraction C, which are indicted in Fig. [3l see 
[35t EH \57\ and references therein. 

• Fermi motion enhancement, C, > 0.8: when the parton carries most of the momentum of the 
nucleon the Fermi motion of the nucleon itself in the nucleus becomes important. 

• EMC effect (proper), 0.3 < ^ < 0.8: the kinematic region of the original discovery, named after 
the experiment, exhibits a suppression Ra{0 < 1 which is usually explained with nuclear binding 
effects. 

• Antishadowing, 0.1 < ^ < 0.3: a region of enhancement of nuclear parton distributions required 
by momentum sum rules. 

• Shadowing, ^ < 0.1: a region of possibly large suppression of parton distributions. It can be 
understood through multiple scattering in the nuclear rest frame, or parton fusion in an infi- 
nite momentum frame. In the deep shadowing (small-^) region this might lead to a color glass 
condensate picture. We refer to [57| for a modern review of models for the shadowing effect. 
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Figure 3: Schematic sketch of the expected behavior of the ratio Ra of nuclear parton 
distributions compared with free nucleon parton distributions. The distinct regions, the 
Fermi motion region, the EMC region and the shadowing and antishadowing regions are 
visible. The ratio Ra is directly reflected in the ratio of single particle spectra in p+A 
collisions to p + p collisions as a function of Pt- 



A parton with 10 GeV/c transverse momentum produced at midrapidities {y = 0) in collisions 
at RHIC energies {^/sNN = 200 GeV) comes from initial parton momentum fractions around C,t = 
2pt/ VSnn = 0.1. Hence it is easy to see that for perturbative calculations at RHIC mostly the 
shadowing and ant i- shadowing regions are of importance. For not too small momentum fractions ^ 
nuclear parton distributions are still in the universal DGLAP regime. They can be measured in deep 
inelastic scattering on nuclei, while the perturbative evolution in the scale fi can be used as a consistency 
check. The parameterizations (which are often parameterizations of the modification Ra for specific 
sets of free nucleon parton distributions) can then be used for hadron-nucleus and nucleus-nucleus 
collisions. DGLAP parameterizations are available from several groups [5H1 E21 EHl ED E21 E21 EH ES] • 
Some, like the EPS08 and EPS09 parameterizations [601 EI], already include some RHIC data in the 
DGLAP fit. This has been done to improve the lack of suitable deep-inelastic scattering data on nuclei. 
Previous deep-inelastic scattering experiments off nuclei cover only large ^ and have very little power 
to constrain the nuclear gluon distribution. This situation leaves us with huge theoretical uncertainties 
on the nuclear gluon distribution below ^ ^ 0.05. Fig. H] shows several fits for modification factors for 
valence quarks, sea quarks and gluons respectively. The spread of possible values for the nuclear gluon 
distribution is truly remarkable. This uncertainty has profound consequences for pQCD predictions at 
LHC energies where the average parton ^ will be much smaller than at RHIC. 

2.2.2 Higher Twist Corrections 

Nuclear corrections to the parton distributions deal with the effects of nuclear binding on the long- 
distance behavior of a process. One can also ask the question whether the hard process between 
partons is affected as well. Indeed it turns out that certain high-twist corrections become important 
in collisions involving nuclei. Corrections beyond leading twist were first characterized in terms of new 
operators beyond parton distributions that appear in the operator production expansion. Twist t was 
defined as the dimension minus the spin of a local operator. The definition can be generalized to apply 
also to situations where an operator product expansion is not available. The leading twist operators in 
parton distributions, qj'^q and F+^F^^, are classified as twist t = 2. We will use higher twist as a simple 
power counting scheme in terms of the large scale Q, such that a twist-t contribution is suppressed by 
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Figure 4: Comparison of different leading order DGLAP parameterizations of the nuclear 
modification i?pb(x, /i) for lead nuclei at /z = 1.3 GeV. The parameterization correspond to 
EKS98 [58J, EPS08 [60j, EPS09 [61], HKN07 [6l], and nDS [ES]. The large theoretical uncer- 
tainties at low momentum fraction x, in particular for the gluons, is clearly demonstrated. 
Figure reprinted from [61] with permission from JHEP. 



1/Q*~^ compared to leading twist. Jaffe's review [66] offers a discussion of both the rigorous and the 
power counting definition of twist. 

Higher twist effects are obviously important if the large scales Q becomes to small (close to non- 
perturbative scales), or if other enhancement effects weaken the power suppression. It was first pointed 
out by Luo, Qiu and Sterman [671 EHl EH] that in large nuclei with mass number A some operators do 
not follow a classification in terms of an expansion in A/Q ^ 1, but rather in the parameter 

AA^/^Q ~ A^L/Q » A/Q . (24) 

Here A is a soft scale (of the order or Aqcd or the constituent quark mass) and L ~ A^^^/A is the 
thickness of the nucleus. L comes into play because in thick nuclear matter multiple hard scattering is 
possible and its probability increases with thickness. Multiple scattering should not modify the total 
cross section very much, but we expect some observables, e.g. transverse momentum spectra, to be 
significantly altered by multiple additional "kicks" that a scattered particle experiences. The Cronin 
effect discussed below is a good example. 

We want to review a simple example, the nuclear Drell-Yan process A + A-^l^ + l^+X [7CTI [7T| 
[72| [73] . At leading order C(a^), and leading twist the virtual photon is produced through a simple 
quark-antiquark annihilation, g + g — > Z+ + /~, see left panel in Fig. [51 The corresponding cross section 
for dilepton pairs of mass Q and (pair) transverse momentum qt is 



da 



dQ'^dq^ 



where 



- anYSiq'r) J2< f [f,/A{^a)f,/B{^,) + fg/A{Uf,/B{^t)] d^, (25) 

J Br, 
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Figure 5: Left panel: Drell-Yan at leading twist and leading order in a^: quark-antiquark 
annihilation. Right panel: A nuclear enhanced twist-4 correction in A + A collisions that 
corresponds to double scattering of the quark from the proton in the nucleus. First the quark 
scatters off a soft gluon and then annihilates with an antiquark. 



essentially is the cross section between partons q and q, Nc = ?>., Cq is the charge of quark q in units of 
e, and = Q'^/i.ibS) is fixed with Ba = Q'^/S. ( |25|) is a straight forward but nevertheless questionable 
result. The g^-spectrum is actually not well defined in the coUinear limit. Indeed the presence of two 
scales qx « Q presents additional problems. The safe way to discuss this result is by using moments 
in g^-space. The lowest moment is the cross section differential with respect to the mass squared 

duAB-^i+i- /"°° da 



dQ^ 7o dQ^dql'"^"^^' ^^^^ 

while the next moment can be used to define the average transverse momentum squared 

M = • (28) 

At leading order and leading twist {qj_) = which is also true for all higher moments. 

The first nuclear enhanced higher twist correction corresponds to double scattering of one of the 
quarks off an additional gluon from the other nucleus, e.g. q + qg ^ + , see right panel in Fig. |5l 
The cross section is 

d(T^AB^i++i- 2\ ^-nassr^ 2 



-5 (gr)^DY-r^2^e, 



dQHql ~ 
1 

[fq/A{Ca)Tgg/B{Cb) + fq/AiCa)Tqg/BiCb) + Tqg/Ai^a) fq/ B^Cb) + Tgg / A^ia) f q/ B^ib)] d^b (29) 

Ba 

Note that there is a derivative on the (5-function. The new matrix elements Tab/ a measure a "soft-hard" 
two-parton density in the nucleus, e.g. 

Tgg/A{0 = I dy^^^e^^^^y^Qiy^ -y^)e{-y^) (A(P)| g(0)7+g(yr) 

X Fa^^iy,)Fa,^y,) |A(P)) (30) 

where P"*" is the large momentum component of a nucleon. Soft-hard in this particular case means that 
the quark or antiquark has a finite momentum fraction ^ while the gluon is very soft. Formally the 
soft-hard matrix elements are limits of more general 2-parton distributions fgg/Ai^qi^g) with 

^..Mfe)= limG4';^(e„ej- (31) 
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Luo, Qiu and Sterman have classified the relevant twist-4 matrix elements that show nuclear en- 
hancement. They all have a probabilistic interpretation as 2-parton densities and they lead to soft-hard 
and hard-hard double scattering on the parton level. The actual nuclear enhancement factor comes 
from unrestricted spatial integrations along the light cone. The coordinates associated with the two 
partons — corresponding to /2 and {y^ + y^)/2 in Eq. (130|) — can be as far apart as the nuclear 
matter extends along the light cone. Parametrically we have 

T,,/^(0 ~ (32) 

where A is a soft scale and A denotes the mass number of nucleus A. 

To arrive at infrared-safe results we again take moments. We note that t = 4 double scattering does 
not make a contribution to the integrated mass spectrum, da^'^'^^ / dQ'^ = 0, or the total cross section. 
However, it leads to non- vanishing transverse momentum, despite the use of coUinear factorization and 
the absence of radiation, 

2\ Airas Y^q^l Ib^ [fq/A{.iq)Tqg/B{.ib) + fq/Ai^a)Tqg / si^b) + ^qg/A (^a) /g/B (^fe) + Tgg/Ai^a) fq/si^b)] d^b 

E,e2/i^ [fq/A{ia)fq/B{ib)+fq/A{U)fqlBm d^b 

(33) 

The t = 4 matrix elements are universal functions that could in principle be measured, but useful 
information is scarce. Most of the time the soft-hard matrix elements are simply modeled using the 
shape of the hard parton distribution 

Tgg/AiO = AM^/V,M(0 (34) 

where A is a parameter with the dimension of energy which parameterizes the strength of the soft gluon 
field. For a symmetric situation with both nuclei being identical this leads to the simple estimate 

(4) « i^2A^ (35) 

Higher twist corrections for t > 4, correspond to multiple scattering beyond double scattering. It 
is possible to identify the operators with maximum nuclear enhancement ~ yl(*-2)/6 ^^^^ they can be 
resummed in certain situations. This is safe to do for Drell-Yan in p+A collisions where the proton can 
be treated at leading twist [731 El]- The resulting effect is a diffusion of qr in transverse momentum 
space. However, generally caution is necessary in nuclear collisions. Although the Drell-Yan process is 
rather simple, with no non-perturbative hadronic structure measured in the final state, factorization still 
breaks down beyond twist t = 4 [2T|. In other words, while nuclear enhanced higher twist corrections 
can be reliably calculated for Drell-Yan in p + A, there are true non-perturbative contributions that 
invalidate this expansion in A -|- A collisions at the level of twist-6. 

Nuclear enhanced higher twist corrections have been considered for several observables, including 
deep-inelastic scattering on nuclei [75], jets and dijets in electron-nucleus collisions [691 US], Drell-Yan 
both at low and high qt [101 [721 El], direct photon production [77j, and photon bremsstrahlung for jets 
[75] . Note that higher twist corrections can appear both as initial and final state interactions. In fact, 
in most cases higher twist corrections could lead to both effects and can not be put in one of those 
two categories. However, those more general cases have not been considered in full detail, and we will 
mostly assume here that higher twist corrections in the initial and final state are independent of each 
other. The most important applications to date for the scope of this article are the Cronin effect (in 
the initial state) which we will discuss next, and medium- modified fragmentation functions (in the final 
state) which will be reviewed in more detail in the next subsection. We conclude by noting that there 
is a patchwork of relevant and useful calculations on the topic of nuclear enhanced higher twist, but a 
lack of comprehensive and systematic studies. 
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2.2.3 Cronin Effect 



The Cronin effect was one of tlie first nuclear modifications discovered in experimental data (TH] . It was 
found that cross sections of hadrons scale with a power of the atomic number A that is larger than 1 
for intermediate transverse momenta Pt ^ 1 GeV/c. This effect was found to not affect the total cross 
sections very much, and to die out like a power law at larger Pt- This is reminiscent of higher twist 
corrections and indeed these results can be interpreted in the framework of higher twist. Intuitively, 
the Cronin effect comes from multiple scatterings of partons on their way to the hard collision. These 
random kicks endow the parton with additional transverse momentum. This leads to a depletion of 
partons with very small (initial) intrinsic transverse momentum and an accumulation of partons at 
intermediate transverse momentum. At even larger values of Pt the additional momentum kicks do not 
play a role and the effect decreases in importance. 

The Cronin effect in its purest form can be studied in the case of dilepton or photon production in 
p+A or A+A collisions. Then it is guaranteed that deviations from the cross sections found in p + p 
are initial state effects. We can simply refer to the discussion from the last subsection where we have 
established that higher twist corrections to the Drell-Yan process that correspond to double scattering 
lead to an increase in the average transverse momentum squared (q^) which is proportional to A^^^, 
and that a resummation of multiple scatterings leads to a Gaussian distribution of qx even at leading 
order in as- 
One can argue that these effects also increase (P^) in hadron production, although it is not always 
clear how to distinguish the effects of initial and final state interactions. It is then quite common to 
refer to less rigorous but phenomenologically successful descriptions of the Cronin effect, see e.g. [80] 
for a review. These models are usually built on the notion of an intrinsic transverse momentum of 
partons in hadrons or nuclei. The concept of intrinsic transverse momentum kj- is not compatible with 
coUinear factorization but has a long history as a phenomenological extension of the former. True 
schemes for A^r-factorization do exist, but only for a handful of select processes, and they are technically 
more complex. Nevertheless many features of the Cronin effect can be described by a model in which 
the average is enhanced in nuclei through 

(kx) pA,AA = {^t)pp "I" ^^T (36) 

where Sk^ scales with the thickness of nucleus A. This is also known as /cT-smearing. 

One possible implementation is through the use of /cr-dependent "parton distributions" 

/aM(e,fcT,/i) = -Lre-^'-'^'-^fa/Aii.li) (37) 

in which {k^) can be fitted to the system size, or calculated from an underlying microscopic model, like 
a Glauber [81] or dipole model [82] ■ Reference [ED] contains a compilation of parameters suitable for 
both RHIC and LHC energies. As a side remark we note that both shadowing and the Cronin effect 
are also natural consequences of gluon saturation and the mechanisms discussed here should smoothly 
transition to their color glass counterparts for very large center of mass energies [83l l8l] . 

2.2.4 Phenomenological Consequences of Initial State Effects 

Initial state effects are considered background effects in heavy ion physics. They are sometimes called 
cold nuclear matter effects, although the two terms are not synonymous. In fact, there are clearly final 
state effects in cold nuclear matter as seen in hadron production m e + A collision by the HERMES 
experiment [85] and successfully described in terms of higher twist corrections [86]. For hadron pro- 
duction and similar process in A + A collisions a factorization between initial and final state effects for 
hadron production is not obvious. It is one of the big assumptions of the hard probes program in heavy 
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Figure 6: Cold nuclear matter effects for pion production in d+An collisions at RHIC. The 
modification factor RdAu for 7r° is shown as a function of transverse momentum P^. Data 
from the PHENIX experiment [8^ is compared to several calculations using different sets 
of parton distributions with or without multiple scattering in the initial state. See text for 
details. Figure reprinted from [87J with permission from Elsevier. 



ion physics that final state effects can be factorized off and modeled separately from hard processes and 
initial state effects. 

One can analyze the effect of initial state interactions in nuclei by looking at hadron production in 
p+A and d+A collisions, and by studying photons and dileptons for p+A, d+A and A+A collisions. 
While these analyses are not yet complete, the picture that starts to emerge is that for Au ions at RHIC 
energies the initial state effects do not change the yield of hadrons for high Pt > 2 GeV/c by more than 
20%. 

Fig. [6] shows calculations by Levai et al. [87] of the nuclear modification factor 

dN'^^/dPr 

- {N^^n)dNPP/dPr ^^^^ 

for pions in deuteron gold collisions compared to experimental data from PHENIX [88]. (Ncoii) is 
the average number of binary nucleon-nucleon collisions expected for the given centrality class. The 
calculations use different sets of nuclear parton distributions [58l|62], with and without smearing through 
multiple scattering, and are also compared to the transport model HIJING [HH]. The modification 
factor is centered around 1 with very moderate deviations, but we can clearly identify how the regions 
of ant i- shadowing and the EMC effect map onto hadron- P^^ at midrapidity for the RHIC top energy 
of 200 GeV. The calculations using HIJING and HKM nuclear parton distributions [62] together with 
multiple scattering are in reasonable agreement with data except for the most central bin. However, 
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the calculation with only EKS modifications to parton distributions is doing equally well. This suggests 
that higher-twist modifications to hard scattering and modifications to parton distributions are not 
easily separated with the available level of data accuracy. 

This result can be confirmed by looking at the same modification factor Raa for direct photons in 
Au+Au collisions. Skipping ahead to Fig. US] we see that Raa for photons is very close to 1 for Pt > 4: 
GeV/c due to the absence of final state effects. The existing, small deviations can be understood with 
the arsenal of initial state effects discussed in this subsection. We conclude that initial state effects in 
A+A collisions seem to be under control at RHIC energies. 

Nevertheless there are some significant gaps in our understanding going forward to LHC We have 
already mentioned our poor knowledge of the nuclear gluon distribution at smaller values of ^. On a 
deeper level, it is a very difficult task to separate corrections to parton distributions (~ logQ^) from 
higher twist corrections to hard processes (~ Q^) with data covering only a limited amount of phase 
space. A separation into a contribution that follows DGLAP evolution and a power suppressed part has 
large uncertainties. The use of a limited sample of P^-spectra from RHIC for nuclear parton distribution 
fits bears the danger of introducing (erroneously) features into nuclear parton distributions which are 
non-universal. The future Electron Ion Collider should be able to improve this situation tremendously. 

2.3 Final State Effects and Energy Loss 

The main goal of the hard probes program in heavy ion collisions is the determination of basic transport 
properties of the quark gluon plasma. The idea that a hot medium formed in nuclear collisions should 
lead to energy loss of partons in the final state, and to a partial quenching of high-Pj- hadrons, was 
proposed many years ago by Bjorken [5U]. In the 1990s it was realized that the most efficient process of 
energy loss is through induced gluon bremsstrahlung. This topic was covered from several angles in the 
1990s in seminal works which estimated the effect on partons in perturbative plasma [91], for partons 
interacting perturbatively with static scattering centers (GW model) [921 [93] and for multiple soft 
scatterings (BDMPS model) [911 [961 EZ] • Energy loss through elastic scattering had been calculated 
and was generally found to be smaller than radiative energy loss for light quarks and gluons. 

In this subsection we describe some of the underlying concepts and the most important modern 
implementations of parton energy loss. We also comment on some more recent developments including 
jet shapes and jet chemistry. We will focus on light quarks and gluons. We would like to point the 
interested reader to the review by Majumder and Van Leeuwen, recently published in this journal [98] . 
for complementary information, and in particular for a detailed derivation of the higher twist energy 
loss formalism and for a discussion of heavy quark energy loss. 

2.3.1 Basic Phenomenology 

Partons of mass m produced in hard QCD processes are typically off-shell, and the virtuality u = 
\/p^ — w? is on average of the same order as the scale Q of the momentum transfer in the hard process. 
The outgoing parton will radiate bremsstrahlung to get back to the mass shell, producing a parton 
shower and eventually a jet cone. This is an example for vacuum bremsstrahlung. Note that this 
picture is consistent with our earlier discussion of hard processes where large angle radiation in the 
final state would be counted as a higher order correction to the hard process while collinear radiation 
is resummed into fragmentation functions. 

A particle that exchanges momentum with a medium will also change its virtuality with each inter- 
action. It too will radiate bremsstrahlung to get back to the mass shell. This increased rate of radiation 
(or "splitting" ) is an effective mechanism to carry away longitudinal momentum, and it acts as a diffu- 
sion mechanism for transverse momentum (directions are relative to the original particle momentum). 
This additional medium-induced bremsstrahlung obviously depends on the density of the medium, or 
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more precisely on the rate at which additional virtuality can be transferred by the medium. This leads 
to the definition of the transport coefficient 



(39) 



? = 



L X 



which measures the average squared transverse momentum transferred to the particle that propagates 
over a distance L, or equivalently the average momentum transfer squared per interaction, /i^, divided 
by the mean free path A of the particle. 

It was realized early on that destructive interference is a key ingredient of these calculations. This 
is a well-known effect in QED, named after Landau, Pomeranchuk and Migdal (LPM) [99| llOOj . Let us 
consider the emission of a gluon g from a quark q that has an initial energy E. If the relative transverse 
momentum between the partons in the final state is k^- and the energy of the gluon is u then the 
formation time 

tf-^ (40) 

estimates when the final quark-gluon pair can be treated as two independent, incoherent particles. If the 
mean free path A of the quark is of the order of the formation time or smaller, radiation is suppressed. 
In that situation the quark scatters coherently from iVcoh ~ ^coh/A scatterers in the medium. For hght, 
relativistic partons the coherence length is given by the formation time /coh = t/- 

Depending on the energy u of the emitted gluon radiation one can qualitatively distinguish three 
domains for induced radiation in a medium of finite length L |95j : 

• The incoherent regime for small u in which the gluon radiation spectrum is independent of the 
length of the medium and the total energy loss AE would be proportional to the length L. 

• The completely coherent regime for large energies u in which the particle scatters coherently off 
the entire medium and the energy loss AE is independent of the length L. 

• The LPM region in between the two extremes in which scatterings off groups of A^^coh ~ ^coh/A 
particles in the medium are coherent, and several or many of such interactions occur. To determine 
the energy loss the differential gluon spectrum per unit length x, udl/dudx ~ Xj^fu^ has to be 
integrated up to the limit cUcr which corresponds to /coh = i-e. the boundary to the completely 
coherent regime. For any given path length x that critical value is, see Eq. (139 p . 



The LPM effect is expected to dominate the behavior of induced gluon radiation in heavy ion collisions. 
The L^-dependence is a characteristic signature of this effect. 

In the following we will discuss several modern implementations of parton energy loss in more detail. 
They all differ in some of the underlying approximations made. 

• The Higher Twist formalism developed by Guo and Wang [86l 1101] . It derives from the notion of 
higher twist corrections for final state partons in e+A collisions. 

• The AMY formalism, based on the work by Arnold, Moore and Yaffe |102l 1103] . It is based on 
hard thermal loop resummation in a perturbative plasma. 




(41) 



This leads to an energy loss rate 

dE 




(42) 



and an energy loss AE ~ — over the entire length of the medium. 
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• The ASW formalism by Armesto, Salgado and Wiedemann |104^ 1105] which resums multiple soft 
gluon emission as in the BDMPS approach in a finite length medium using Poisson statistics. 

• The GLV approach by Gyulassy, Levai and Vitev |106[ I107[ 1108] which considers hard scatterings 
off static scattering centers in an opacity expansion. 



2.3.2 The Higher Twist Formalism 

The systematic discussion of final state interactions of hard scattered partons in a nuclear medium 
is dominated by several big questions. One of the most fundamental ones is whether the final state 
interactions be factorized off (a) the hard process, (b) the initial-state effects in nuclei, and (c) the 
fragmentation into hadrons? There are ways to treat problem (c), or it can be circumvented by looking 
at jets instead of hadrons, which is experimentally difficult at RHIC, but will be routinely done at the 
LHC. Most of the QCD-inspired energy loss models that we discuss here assume such a factorization. 
The Higher Twist (HT) formalism eventually has to make the same assumption, but it takes guidance 
from a process in which such a factorization can actually be tested: semi-inclusive hadron production 
in deep inelastic scattering e+A off nuclei. 

Guo and Wang were the first to write down a set of expanded evolution equations for medium- 
modified fragmentation functions in e+A collisions [86| 1101] . They base their computation on the 
pioneering work of Qiu and Sterman on nuclear enhanced higher twist corrections discussed earlier. 
Semi-inclusive hadron production, e + A^e + H + Xis usually discussed by factorizing the cross 
section o" as a function of hadron momentum l'^ = {E^, in) and final lepton momentum = {E2, P2) 
into a QED part called the leptonic tensor and a QCD part called the hadronic tensor 

This factorization is accurate to leading order in the electromagnetic coupling aem- The leptonic tensor 
is 

L^, = 2p1pl^ + 2p'ip'i - 2pi ■ p2 g^u , (44) 

and = —q^ measures the virtuality of the photon with momentum = P2 —Pi- We have labeled the 
initial momentum of the lepton as p'^ and as usual we call the average momentum of a nucleon in the 
nucleus with a large light cone momentum fraction . It is common to choose the frame such that 
the photon momentum has a large — component and no transverse components, q^ = {—Q'^/2q~, q~, 0) 
in light cone notation. 

The hadronic tensor measures the electromagnetic current (to which the photon couples) both in 
the amplitude and complex conjugated amplitude between initial nuclear states and the final hadronic 
states, W'^'^ ~ ^■^(Alj^lif, XIj'^IA). After integrating the transverse degrees of freedom we can 

write it as a function of only the longitudinal momentum fraction of the hadron with respect to the 
photon momentum, zh = IJi/q' ■ Leading-twist {t = 2) coUinear factorized QCD tells us that 



E / «M(e,Q) —H^^{^,p,q,z,)D,/H (^,Q 



(45) 



where Hj^l^ describes the hard scattering of parton a off the virtual photon, producing a parton b and 
maybe more unobserved final state particles x, a + 7* — > b+x. While parton a has a momentum fraction 
^ with respect to P^, parton b has a momentum fraction Zb = lb /^^ with respect to the photon. There 
is a clear separation between the initial and final state long-distance processes described by the parton 
distribution / and the fragmentation function D respectively. The leading order diagram for this process 
is shown in the left panel of Fig. [7] 
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Figure 7: Left panel: Semi-exclusive deep-inelastic scattering at leading order. Right panel: 
example for a real radiative correction that leads to the evolution equations for the fragmen- 
tation functions. 




Figure 8: Two examples of twist-4 diagrams that contribute to the evolution equations of 
medium-modified fragmentation functions. Propagators with poles that make the gluon from 
the nucleus soft (momentum ~ ^^i) are indicated by a grey circle, propagators whose poles 
make the gluon harder (momentum ~ are shown with a solid black circle. There are 
many more diagrams including those with different final state cuts and interference diagrams 
in which the gluon from the nucleus couples to different particles in amplitude and complex 
conjugated amplitude. 



Collinear radiation off the final state parton b leads to leading-twist evolution equations for the 
fragmentation functions. The diagram in the right panel of Fig. [7] leads to a correction to Eq. f H5|) 
which can be written as 

with the familiar splitting functions Pt^c to a third parton c which undergoes fragmentation. The 
leading collinear term can be resummed and leads to evolution equations which are completely analogous 
to the DGLAP equations ([7]) for parton distribution^. 

Guo and Wang showed that the hadronic tensor at the level of nuclear enhanced twist-4 receives 
contributions from diagrams like the ones shown in Fig. [H] where parton b or c, or the radiated parton 
could scatter off an additional medium parton d. They computed the result from those diagrams and 



^The equations in [5511101] often omit the integral over Zf, since the hard parton scattering tensor H contains a 5{zi, — l) 
function at leading order in only, which has been canceled with the integral. 
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(47) 



For simplicity of notation we have assumed here that the second parton d is a gluon that does not 
change the identity of the parton it couples to. The other cases can be treated accordingly. The twist-4 
matrix elements are similar to the soft-hard matrix elements introduced in Eq. fl30|) . If parton a is a 
quark we have 



Tgg/A{^,^L) 



X (A(P)I g(0)7+g(z/r)i^r (l/3~)i^a/(l/r) • (4J 



The expression also contains a derived matrix element 



dz 
l-z 



(49) 



which remains after the virtual corrections have canceled the singularity at ^ = 1. 

We note that the normalization of the matrix elements, that is following Guo and Wang here, differs 
by a factor 2tt from Eq. fl30|) . We have also gone beyond the soft-hard matrix elements by allowing 
parton d to have a non- vanishing momentum fraction 



Jjrp 



2P+q-z{l - z) 



(50) 



with z = Zc/ Zh. It is the structure (1 — exp(. . .))(! — exp(. . .)) in PHI) that exhibits interference and 
will eventually lead to the LPM effect in the Higher Twist formalism. One can easily see how this 
interference emerges in the calculation. The momentum of parton d is fixed by the poles of partons b/c 
and there are two kinematic possibilities shown in Fig. [HI Either the momentum of parton d is very soft 
with a momentum fraction C,d ~ k^/Q"^ which vanishes when its intrinsic transverse momentum is set 
to zero. The phase exp{—iC,DP^y~) then reduces to unity. The other pole sets the momentum fraction 
to xl, and interestingly the amplitudes for both poles exhibit a relative minus sign. This interference is 
common when higher twist corrections are considered together with radiative corrections. A discussion 
of this soft-hard interference in the context of Drell-Yan can be found in [70[ [71], [721 1109] . 

The collinear radiation corrections at twist-4 can be resummed in modified evolution equations for 
new fragmentation functions D in cold nuclei just as in the twist-2 (DGLAP) case. The resulting set of 
equations takes exactly the same form as in Eq. ([7]), with D ^ D and new, medium-dependent splitting 
functions 

Pabiz) = Pabiz) + APabiz) (51) 

where Pab is the usual vacuum splitting function and APab is the twist-4 correction. For the case of 
q q + g splitting induced by a gluon we have 



27TasCA 1 



Jjrp 



l+Z^\ Tgg/A{^,i 



1- z 



f,/AiO 



+ S{z-1 



T) 



(52) 



The other modified splitting functions are discussed in Ref. [101] for scattering off gluons and in Ref. 
[110] for scattering off quarks. Even though this result is technically correct and very useful we can also 
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see its limitations. The splitting functions, and as a consequence the modified fragmentation functions 
D, are no longer universal as they depend on the underlying process through the matrix elements T. In 
fact they are no longer just functions of a single momentum fraction z but also depend in a non-trivial 
way on ^. On the other hand, the breaking of universality encodes the medium effects that we are after. 

The new medium-modified fragmentation functions allow us to write hadron production in e+A 
collisions in a very simple way analogous to Eq. P5|) as 

^ = E / dUa/A{i. Q) C — g, z,)b,/H (—, q] . (53) 

dZH ^ J Jz„ Zb \Zb / 

The dependence of the D on other quantities is usually suppressed in the notation. We have now come 
to a point where one has to introduce a certain amount of modeling since a rigorous solution would 
include a simultaneous fit of the Dab/A and the Tab/ a in the same environment (because of the loss of 
universality) with the evolution equations as constraints. This is too complex a task given the available 
data. 

The twist-4 matrix elements are modeled similar to the less general soft-hard matrix elements from 
Eq. ( ETj) . It seems safe to assume that Tab/ a can be factorized into a product of two parton distributions 
for partons a and h resp. Guo and Wang model the interference effect by introducing a massless 
parameter for the radius Ra of the nucleus, xa = ^/{MRa), where M is the mass of a nucleon. They 
suggest 

Tab/Ait a) ^ — (l - e-^l/^-) [fa/AiO Ub/A{U + fa/A^i + ^l) >^b] (54) 

Xa ^ ^ 

where C is a normalization constant and Hb ~ linix^o xfb/A{x) formally is the value of the parton density 
for parton h when it is very soft, i.e. with momentum fraction of order xd- Obviously Kg is proportional 
to the gluon density introduced in Eq. ( IMll . Note that the dependence on the size of the system is 
hidden in xa and the leading size dependence is ~ A^/^. On the other hand the formation time of the 
radiation is ~ 1/{xlP^) and the factor 1 — e'^^I^A leads to LPM suppression unless ^ xa- In Ref. 
|112j the authors suggest an even simpler model that drops the second term 



Tab/A{^, a) ^ ^ (l - e'^l/^'A fa/AiO ■ (55) 

Xa ^ ^ 



Even if the set of twist-4 matrix elements were perfectly known there is still the task to extract the 
medium-modified fragmentation functions from the set of evolution equations. Guo and Wang originally 
suggested the first iteration 



- as /"^^ dl^ dy ^-^ f z 

Da/H{z,fi) = Da/H{z,fi) + TT -JT — ^Pa~,b{y) Db/ H - 



(56) 



as an approximate solution. Note that the APa-^b carry all the information about the medium through 
Eq. (152!) and its counterparts. Deng and Wang recently showed how to solve modified evolution equations 
numerically (111] . 

Using the first iteration the energy loss of a parton can be calculated as the shift in the momentum 
fraction due to the term APa^b in the equation above 

(Ay) = l-Y r ^ ['dyJ2 ^Pa^biy) ■ (5T) 

JO T J z ^ 

For quarks one finds [86l 1112] 
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which is proportional to the nuclear size squared, (Ay) oc A^/^, as expected for the LPM regime. 
xb = Q^/{'^P^Q~) is the Bjorken variable in deep inelastic scattering. In Ref. [112] the authors can 
explain the observed suppression of semi-inclusive hadron production in the HERMES experiment |113] 
very well using the medium modified fragmentation functions described above. One can go one step 
further and try to interpret the medium modification as a simple rescaling of the vacuum fragmentation 
functions. One uses an ansatz |112l I114j 



where A2 is the typical energy loss for parton a. This formula can only be a satisfying approximation 
in a limited region with z + /\z < 1. At this level the medium-modified fragmentation functions have 
been cast in a very general form and one can try to apply the general concepts to systems other than 
e+A. In particular one can extract the stopping power (AE/L) EAz/L for a particle with initial 
energy E from data in heavy ion collisions. Using the techniques here (and including the dilution of the 
medium through the longitudinal expansion in nuclear collisions) E. Wang and X.N. Wang concluded 
in 2002 that the differential energy loss extracted from data in Au+Au collisions at RHIC was about 
15 times larger than the one for cold nuclear matter measured at HERMES [112] . 

Beyond the first iteration, Deng and Wang have systematically studied the evolution of gluon induced 
parton-to-parton fragmentation functions (starting from Da/a = ^{z — 1) and Da/b = 0, b ^ a initial 
conditions at a low scale). They also investigate the energy and medium thickness dependence and 
apply their techniques to pion fragmentation with parameterized vacuum fragmentation functions as 
input. Modified evolution suppresses the fragmentation functions at intermediate and large values of z 
as expected [111] . 

In recent years progress was made on the HT formalism (in its original meaning for deep-inelastic 
scattering) by considering medium modifications for double fragmentation [ 115] , elastic energy loss [116] 
and through successful resummation of multiple scatterings per photon or gluon emission [THl 1117] . 

2.3.3 The AMY Formalism 

The particular merit of the formalism widely know as AMY, after Arnold, Moore and Yaffe, is its 
complete internal consistency in a very high temperature regime. The basic picture is the following. 
Partons propagating through a plasma with temperature T, themselves having momenta of order T 
or larger but small virtualities, interact perturbatively with quarks and gluons in the plasma with 
thermal masses ~ gT. They pick up transverse momentum of the same order gT, and then radiate 
gluons (for quarks) or split into quark- ant iquark pairs or two gluons (for gluons), again at typical 
transverse momentum scales gT. This leads to formation times that are rather long, of the order 
T/{gTY ~ ^/{q'^T). The shortcomings are the requirements of small initial off-shellness of the parton, 
an unlikely condition for a parton emerging from a hard process, and the condition of very small coupling 
(7 ^ 1 to justify thermal perturbation theory, which requires very large temperatures T ^ T^. The 
exact temperature from which on such calculations start to be reliable is a matter of debate. One 
must also assume that the temperature does not change during the formation time of radiation which 
is questionable for rapidly evolving fireballs. Nevertheless, the rigor of the formalism has made AMY 
an appealing choice in the canon of energy loss calculations. 

The AMY formalism grew out of a computation of the complete leading order, hard thermal loop 
(HTL) resummed perturbative photon and gluon emission rates. Arnold, Moore and Yaffe introduced, 
for the first time, the correct treatment of coUinear emission in a finite temperature medium, which 
must take into account the LPM effect due to the long formation times ~ l/{g'^T) [511 [521 [102]. The 
applications of these results to photon production in quark gluon plasma had been mentioned in a 
previous section. The results for gluon radiation off a parton with typical momentum T ^ gT in the 




(59) 



25 



medium can be used to calculate its rate of energy loss. For quarks of momentum p radiating gluons of 
momentum k the rate is |1U2[ llU3j 



dkdt 



(P^P-k) = ^ns{k)nAp-k)'-^ll-^ J A^2h.ReF(h,p,A:) (60) 



where hb and np are the boson and fermion occupation factors for the gluon and the quark in the 
final state, x = k/p is the momentum fraction of the gluon, and h = p x k is a useful measure for the 
non-coUinearity of the final state. F oc h is a function defined by the integral equation 

2h = i6E{h,p,k)F{h) + ^ [ d'q^Viq^) x [(2Cr - Ca) (F(h) - F(h - A;qx)) 

2iT J 

+ Ca (F(h) - F(h + pq^)) + Ca (F(h) - F(h - - A;)q^))] . (61) 

We have dropped the second and third argument in F{h, p,k) for brevity, in all cases above they are 
equal to the initial momentum p and the radiated momentum k resp. It is the function F that encodes 
important properties of the medium via the potentials 

as a function of momentum transfer q^ and temperature T. 

mh,P,k)= ,,+^+ (63) 

^ ^' ^ 2pk{p-k) 2k 2{p-k) 2p ^ ' 

is the energy difference between final and initial state. The masses m are m£,/ \/2 for gluons and gT / -\/3 
for quarks in the respective channels with momenta p, k and p — k. 

The rates for other processes can be obtained from Eq. (160 p by replacing the splitting functions, 
adjusting the Bose or Fermi factors, and putting the correct color factor Cr. The missing splitting 
functions needed are 



^/-:57T ^ ioig^qq, (64) 



+ (1 — x^"^ 
x'^{l — x^ 
l + x^ + (l-x)^ 
x^(l — x^'^ 



for g^gg. (65) 



The rates for different processes as a function of time t can be implemented in coupled rate equations 
for quark, antiquark, and gluon momentum distributions fq, fg and fg resp. 



dt 
dt 



oo 

dk 

— oo 
oo 

dk 

-oo 



dTq^ggip + k,p) ^^^^^ ^^^^ ^ _ ^WP.fe) ^^(^) 



(66) 



dkdt '•"'''^ ' •'"''^ " dkdt 

(67) 



_^9r,^>^^+M) ^^^^ + ^) - ^^7km^^ ^'^P + - 



The equation for the antiquark distribution is analogous to the equation for the quark distribution. 
Note that the emitted momentum k can be positive or negative, which means that a parton can in 
principle also acquire momentum from the medium. Of course, a parton with momentum much larger 
than typical thermal momenta will still lose momentum on average. Final quark and gluon spectra can 
be subjected to vacuum fragmentation to compute final hadron spectra. 
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2.3.4 The GLV Formalism 



The GLV energy loss model by Gyulassy, Levai and Vitev |106lll07lll08j is based on the earlier Gyulassy- 
Wang model. It describes the medium as an ensemble of static scattering centers with Yukawa potentials 
exhibiting a screening mass fi, which also set the typical scale of transverse momentum transfer. The 
scattering amplitude is then expanded in terms of the opacity L/X where L is the length of the medium 
and A is the mean free path. The leading, zeroth order term, for a parton with energy E corresponds 
to vacuum radiation with a spectrum 

as a function of the gluon longitudinal momentum fraction x and transverse momentum k. Cr is the 
color Casimir factor in the appropriate representation, Cp = 4/3 for quarks and Ca = 3 for gluons. 

At the next order in opacity one considers the interference of vacuum radiation and a single medium- 
induced radiation with momentum transfer q |108] 

{2-2x + x)-E- / dq ,^^ ^ ..2^2^a^2T,2 , ^ur2 - (69) 



dxdk'^ 2n ^ ' \ k"^ j ^ Ti{q^ + ^I'^f IQx'^E'^ + c^fL'^ 

The integrals can be evaluated analytically away from the extreme cases x — )■ and x — )■ 1. Then 
maximally loose kinematic constraints 0<A;<oo,0<g<oo can be assumed and the total energy 
loss at first order in opacity is 

AEW = ^^L2ln-. (70) 
4 A /i 

It exhibits the characteristic L^- dependence. In contrast the zeroth order gluon spectrum fl55|) will lead 
to a "vacuum quenching" 

AEW = ^^Eln- (71) 
37r yU 

where in this case is chosen to play the role of a lower cutoff for the transverse momentum k. 

Higher orders in the opacity (L/A)" can be treated numerically |108] . Since the number of emitted 
gluons is finite, fluctuations around a given mean value are important. They can be taken into account 
through Poisson statistics |118] . leading to a probability distribution P(e) for the fractional energy loss 
e = AE/E. The probabilities for emission of n gluons are 

p-Ng r dT f \ 

where Ng is the average number of gluons and dl/dx is the gluon energy spectrum to the desired order 
in opacity, e.g. derived from Eq. ( I69l) to first order. The total probability distribution is 



Pie) = J2P^ie). (73) 

n=l 

It can be used to define a medium-modified fragmentation function 

^ (74) 



analogous to Eq. ( |59|) . Recall that we tacitly assume that parton energy loss and the actual hadroniza- 
tion of the parton are factorizable, and hadronization itself happens outside of the medium. There have 
also been attempts to include elastic scattering consistently in the GLV approach |119[ 1120] . 
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2.3.5 The ASW Formalism 



The energy loss model due to Armesto, Salgado and Wiedemann [104^ 1105] assumes a Poisson-like 
distribution of gluon emissions as described in Eqs. fl72]) and fITH]) . The proponents of the model present 
a resummed version of the so-called quenching weight P. They add the explicit possibility Pq of zero 
gluon emissions, i.e. a finite probability that a parton escapes unquenched. For their main result they 
use the gluon radiation spectrum dl/dx from the BDMPS approach assuming finite propagation length 
L [nH ESI ESI EZl 1121] . They also present results for a resummation of gluon spectra in a GLV-like 
opacity expansion [122], but we will not discuss that latter option in detail. 

In the case of BDMPS soft scattering they introduce a characteristic gluon frequency 

= ^qL" (75) 

and a dimensionless quantity 

R = UcL. (76) 

Note that Uc is close to the definition in Eq. (14T!) . R is introduced to enforce the kinematic constraint 
k± < u. From its definition we can infer that i? — ?■ oo corresponds to an infinitely large medium if Uc 
is finite. It is also the limit in which the previous BDMPS result is recovered. 

One can perform a numerical resummation of the Poisson sum P from Eq. ( 1731) for > 0. The total 
probability is then written as a sum 

P{AE) = poioo,, R)6{AE) + p{AE, oo,, R) (77) 

which contains a discrete probability po for zero energy loss and the resummed probability for finite 
energy loss (n > 0). The unphysical case AE > E can be dealt with by either renormalizing the 
total probability to unity ( "reweighted" ) or by introducing a second 5-function at AE = E which 
accumulates the probability for total loss of the jet ( "non- reweighted" ) [1211 1123] . The uncertainty in 
the treatment of the case AE > leads to a rather large uncertainty in describing the data, see e.g. 
Fig. [22l The authors of the ASW model provide a Fortran code which computes both the discrete and 
the continuous part of the quenching weight as functions of AE/u, and R [124] . As in the GLV case 
the quenching weights can be used to define modified fragmentation functions using Eq. fl7^ . although 
also in the ASW case the true, non-perturbative fragmentation process itself is considered independent 
and not affected by the medium. 

For applications in heavy ion collisions we need a way to go beyond the simple assumptions of a 
homogeneous medium in which q is constant along the propagation path of the parton. In fact in a 
realistic fireball q = (x, t) where x is the position in the fireball and t the time. The extraction of 
q = (x, t), or at least a spatially and time-averaged version (g) of it, is actually the goal of the hard 
probes program. One can deduce cUc and R from the two lowest moments of q integrated over the path 
of a jet particle, 

poo 

4= / g(x(t),t)|x(t)-xorrft. (78) 







Here x(t) is the trajectory of the particle originating from a point xq at t = (the integrals can also be 
shifted by a finite formation time). Then we have 

2/2 

ujc = h, R=-r- (79) 

-'0 

Due to the sheer impossibility to extract detailed space-time information on q it has become stan- 
dard to assume that q is locally proportional to a quantity whose distribution and time-evolution is 
approximately known. A popular choice is the 3/4th power of the local energy density e 

g(x,t) = 2ire^/^(x,t), (80) 
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or the entropy density s. This is the parametric dependence expected from dimensionahty arguments 
for a fully thermalized quark gluon plasma |125] . K is then treated as a fit parameter and we expect 
it to be close to unity for a weakly coupled plasma and larger for a strongly coupled system. Other 
choices for modeling the shape of q found in the literature are the temperature T(x, t) of an equilibrated 
plasma, and the density of the number of participant nucleons or number of binary nucleon-nucleon 
collisions, which are then usually treated as time-independent medium distributions. 

2.3.6 Final State Effects: Other Developments 

We want to end the discussion of final state effects in nuclear collisions by briefly touching upon two 
special topics. We have mostly focused our attention on the production of hadrons since this is the 
dominant mode of measurement at RHIC. At LHC, the calorimetric measurements of jets will become 
much more important as their energy grows much above the background event. The advent of fast and 
reliable jet algorithms for high-multiplicity environments |126l 11271 1128] has added to the excitement. 

One possible way of modeling jets is through advanced Monte Carlo simulations of medium induced 
gluon radiation that does not just focus on the leading parton but tracks the evolution of the entire 
parton shower. Monte Carlo jet quenching modules like PYQUEN [129], Q-PYTHIA pLSOj, JEWEL 
[131] , YaJEM |132] and MARTINI |133] have made progress in that direction. Another approach is the 
study of jet shapes in heavy ion environments that track the flow of energy through cones of given radius 
[134] . These studies will be much more flexible and comprehensive than leading hadron studies as they 
give much better answers to the question of "energy loss" which is, of course, rather a redistribution of 
energy than a real loss. 

Another interesting field that has emerged in recent years is the study of hadron chemistry in jets. A 
heavy ion environment not only redistributes energies of energetic particles, it can also lead to significant 
changes in the relative abundances of hadrons. This can either happen through profound changes in 
the way hadronization works in high multiplicity environments |135] . see also the discussion on quark 
recombination in the next section, or through the exchange of particles with the quark gluon plasma 
which leads to a phenomenon termed jet conversion |136l 1137] 11381 [T39] . Jet conversions would increase 
the number of protons and kaons relative to pions in nuclear collisions vs p + p collisions. In the former 
case constantly occurring conversions between quarks and gluons wash out the different color factors 
for their respective suppression, leading to equal suppression of particles from light quark and gluon 
fragmentation. In the later case the small sample of high-pT- strange quarks is tremendously enhanced 
relative to up and down quarks when quenched in a chemically equilibrated quark gluon plasma. In 
Ref. [137] the authors predicted a factor 2 increase in the Raa of kaons vs pions which seems to bee 
seen in preliminary STAR data [140] . 

Conversions have been particularly well understood in the case of photons [14H I142[ 1143] and 
dileptons [144[ 1145] . Both induced photon bremsstrahlung [SH [HS] and elastic annihilation and Compton 
scattering with the medium {q + q ^ J + g and q + g — )■ 7 + <? resp.) lead to photon yields that 
are comparable with other sources (thermal, hard direct, vacuum bremsstrahlung) at intermediate 
transverse momenta of a few GeV/c. Both the evolution equations of the Higher Twist formalism and 
the rate equation of AMY can accommodate more channels and "flavor" changing processes in a straight 
forward manner to study these effects. 

2.4 The Perturbative Approach: Critique and Challenges 

Despite the large amount of effort put into the development of a perturbative description of hadron 
production in heavy ion collisions, there are uncertainties remaining about the exact nature of jet- 
medium interactions in the kinematic and temperature regimes important at RHIC. We will discuss in 
more detail in Sec. H] below that the four approaches described here generally fare well in describing 
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Table 1: pQCD-based energy loss models: This table summarizes some of the 
key assumptions of the four perturbative calculations discussed here. The mod- 
els differ with respect to the medium (thermalized, perturbative), kinematics, 
scales {E = energy of the parton, = transverse momentum of the emitted 
gluon, fi = typical transverse momentum picked up from the medium, T = 
temperature, A = typical momentum scale of a (non-thermalized) medium, x 
= typical momentum fraction of the emitted gluon), and the treatment of the 
resummation. 



Model 


Assumptions about the Medium 


Scales 


Resummation 


GLV 


SLdLic scdLLeiiug cenieis i UKiiwa j , opaciLy 
expansion 


^ > /ct ~ /i, 2; < 1 


Poisson 


ASW 


static scattering centers, multiple soft scat- 
tering (harmonic oscillator approximation) 


E ^ kT ^ X <^ 1 


Poisson 


HT 


arbitrary matrix element at scale A (thermal- 
ized or non-thermalized medium) 


-E > A;^ > A ~ /i 


DGLAP 


AMY 


perturbative, thermal, g « 1 (asymptoti- 
cally large T) 


E>T:^ gT ^ fi 


Fokker-Planck 



RHIC data, but they can reach very different quantitative conclusions about the quenching strength q. 
This should not come as a big surprise since the approaches differ in some of their basic assumptions, 
and there are large uncertainties in modeling hard probes beyond the calculation of an energy loss rate 
for a quark or gluon. 

Currently the big picture can be summarized as follows: perturbative calculations under various 
assumptions are compatible with RHIC data, but the constraints are insufficient to rule out any of 
the models. The experimental constraints are also insufficient to completely exclude non-perturbative 
mechanisms of jet quenching. Calculations using the AdS/CFT correspondence to model strongly 
interacting QCD [1461 \1A7\ 1148] can describe the same basic phenomenology. Most likely this challenge 
to perturbative QCD can only be answered at LHC. The extrapolation of jet quenching to larger jet 
energies is significantly different in strong coupling and perturbative scenarios |149j . It is also possible 
to imagine a small regime of strong non-perturbative quenching around Tc together with perturbative 
quenching at higher temperatures. Such mixed scenarios might be hard to distinguish experimentally. 
One such picture was recently explored by Liao and Shuryak |150] . They found that a "shell" -like 
quenching profile in which quenching is enhanced around Tc can give better simultaneous fits to single 
hadron suppression and elliptic flow. 

Additional uncertainties come from a variety of issues regarding the details of the phenomenological 
modeling: 

• The initial state: Jets will interact with their environment before a quark gluon plasma is fully 
formed. The time dependence of q during the first fm/c of the collisions is highly uncertain, in 
particular if initial interactions are dominated by coherent gluon fields fl9\ [20] . Some calcula- 
tions set an ad-hoc start time or use different extrapolations to small times. One estimate of 
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uncertainties can be found in [151] and is discussed in more detail in Fig. HO] below. 

• Fireball evolution: Wildly different fireball parameterizations used in jet quenching calculations 
can be found in the literature up to this day. The correct longitudinal and transverse expansion 
with the correct cooling rate have to be taken into account. Recently comparisons of differ- 
ent calculations using the same underlying fireball calculated from hydrodynamics have become 
available |152j . However, as we will discuss in detail in the Sec. [3l uncertainties remain in hydro- 
dynamic calculations as well which are transferred to hard probes when hydrodynamics is used 

background. 

• The hadronic phase: The models referenced here deal with parton energy loss in a partonic 
medium. Clearly a jet will also interact with a surrounding hadronic medium. Some models 
could in principle deal with this situation (e.g. the higher twist approach only needs the jet to 
be dominated by a sufficiently high energy parton), others have to fail (like the AMY approach). 
But none of them addresses the question of fully formed high energy hadrons in a jet interacting 
with a hadronic medium. Shower simulations with full space-time evolution might be helpful to 
constrain at least the size of the problem from the partonic side. 

• Event-by- event fluctuations: Not much is known from experiment about spatial fluctuations in 
the fireball, but clearly we should expect a fireball to exhibit a certain degree of inhomogeneity 
as suggested by many models of initial nucleus-nucleus interactions. Compared to quenching in 
a smooth, average fireball, event-by-event quenching can lead to considerably different results for 
hadron suppression and elliptic flow |153j . 

• Back-reaction from the medium: While the medium modifies jets, jets on the other hand modify 
the surrounding medium by transferring energy and momentum. The heating of the medium can 
be considerable |154] . and a variety of shock phenomena can occur 1156^ 1157] . Clearly, a 
comprehensive approach will consider both the jet and the medium as variable and would not 
fix one or the other as a background or boundary condition. In particular, if part of the energy 
(or maybe most of the energy for some jets) thermalizes, this most likely proceeds through non- 
perturbative channels which are not included in either of the models discussed here. 

There are no systematic studies of all uncertainties together. A pessimistic estimate of their compounded 
effect would be a factor 3-5 uncertainty in the extraction of q from RHIC data. 

Let us finally revisit the four approaches to calculate leading particle energy loss discussed here. 
Why do they lead to similar qualitative, but sometimes quite different quantitative results? We have 
discussed the underlying assumptions of each model in its respective section. We summarize them 
once more in Tab. [T]in terms of the different ways the medium is modeled, the hierarchy of scales and 
the way resummation is handled. Two key points are the different assumptions about the transverse 
momentum kr of the emitted gluon vs the transverse momentum /i picked up from the medium, and the 
philosophy of multiple soft emissions vs an opacity expansion (with single, somewhat larger transverse 
kicks). The full solution of the problem is quite complex, even in a fully perturbative medium. We refer 
the interested reader to the recent assessment by Arnold |158] and the review by Majumder and Van 
Leeuwen [98]. It might be a while until more comprehensive calculations are available, but efforts in 
this direction are under way, e.g. within the TECHQM and JET collaborations. Simplistic assumption 
have to be improved and narrow kinematic regimes have to be widened. 

3 Success of Hydrodynamic Models at RHIC 

Since Landau |159] and Bjorken jl60] proposed the idea to apply hydrodynamics to the production of 
particles in high energy collisions, it has evolved into one of the most useful phenomenological models for 
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our understanding of high energy heavy ion colhsions. Starting with Landau's hydrodynamic description 
and Bjorken's scahng solution, various kinds of implementations of hydrodynamics were developed to 
understand experimental data from the Alternating Gradient Synchrotron (AGS) and and the Super 
Proton Synchrotron (SPS). Hydrodynamics and hydro- inspired models (e.g. the Blast- Wave Model 
[161] ) give us reasonable explanations for a large amount of experimental data: single particle spectra 
(with respect to transverse momentum Pt and pseudorapidity 77), two particle correlations, collective 
flow, and electromagnetic probes. However they did not appear to work for some aspects of anisotropic 
flow, e.g. directed flow vi and anisotropic flow t>2. For example, the rapidity dependence of directed 
flow of charged pions is different from that of protons [1621 1163] . The charged pion vi decreases with 
rapidity, but the proton vi increases (see Figs. 18 and 19 in Ref. [163] ). This difference between pions 
and protons is difficult to understand from collectivity arguments which is one of the crucial features of 
hydrodynamics. To explain this interesting behavior in detail, hydrodynamics may be too simplistic and 
additional effects may play a role. Therefore transport models in which more complicated dynamics of 
the underlying theory are included fare better for directed flow. Another hint comes from the fact that 
hydrodynamic models routinely overestimate the size of elliptic flow at SPS energies [163] . This suggests 
that the system does not completely equilibrate. Hence the validity of hydrodynamics at SPS energies 
is not very clear. This led to the pre-RHIC view that hydrodynamic models were rather simple-minded 
phenomenological tools. 

All of this changed dramatically after the flrst experimental results from RHIC came out in 2000. 
Hydrodynamic models could naturally explain the unexpectedly large elliptic flow at RHIC compared 
to that at SPS |164[ 1165] . The success of hydrodynamics at RHIC, together with observations of 
jet quenching in the medium many times larger than that in cold nuclear matter, indications for the 
existence of a Color Glass Condensate and the success story of recombination models, are cited as 
the main results from the RHIC program, as evidence for production of a quark gluon plasma with 
deconflned partons, and for the conjecture that this QGP is strongly coupled (a "sQGP") and in fact 
may be the most perfect liquid ever seen [3l 1166] . 

In the wake of this success hydrodynamics has become the most promising tool to describe most of 
the expansion and cooling process of the bulk of the matter created in heavy ion collisions at RHIC. 
However, at the same time limitations of hydrodynamic models, old ones and novel ones, were found: 
e.g. the failure to describe experimental data on two-particle correlations, or elliptic flow as a function 
of pseudorapidity. This suggests that the assumption of perfect fluidity is not valid for the entire flreball 
at RHIC. At least in the hadronic phase and in the flnal state viscous effects can not be neglected [167 . 
Since the perfect fluidity hypothesis is tested through precision analyses of elliptic flow |168] one needs 
to overcome this obstacles with hybrid models that couple hydrodynamics to hadron-based transport 
models and through the development of viscous hydrodynamic codes. 

Both of those avenues have undergone remarkable developments in recent years. In particular our 
understanding of second order viscous, relativistic hydrodynamics has increased tremendously over the 
past few years. While the extraction of the equation of state for quark gluon plasma and the phase 
transition used to be at the forefront of hydrodynamic modeling, the extraction of transport coefficients, 
like shear and bulk viscosity, of the hot and dense matter created at RHIC has been added to the main 
goals. In this section we will review the basic principles and numerical concepts of hydrodynamics 
together with some recent progress. We will then discuss applications of hydrodynamic models at 
RHIC. 
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3.1 Basics of Relativistic Hydrodynamics 

3.1.1 The Framework of Ideal Hydrodynamics 

Let us start with the ideal relativistic hydrodynamic equations of motion. In ideal hydrodynamics 
thermalization is perfect and there is a well-defined local rest frame for each fiuid cell. In that local 
rest frame, the energy-momentum tensor of a volume element of the fiuid (where Pascal's law works) is 
given by 

f^'^(x) 

where e{x) and p{x) are energy density and pressure, respectively, as functions of the space-time point 
x'^ of the fiuid cell. Introducing the four-velocity u^{x) of each fiuid cell in the lab frame we can boost 
the local energy-momentum tensor into the laboratory frame which gives us the well-known result 

Tf^'ix) = [e{x) + p{x)]u^'{x)u''{x) - p{x)g>"' . (82) 

The motion of the fiuid is simply described by the equations of energy and momentum conservation, 

a^T'^'^(x) = 0, (83) 

from which the entropy current conservation, 

d^{su^) = (84) 

is derived. Other conserved currents j^{x) besides energy and momentum — most importantly the 
baryon current J^(a;) — can be taken into account by imposing the conservation laws 

d^3^{x) = Q. (85) 

In ideal hydrodynamics each of these conserved currents can be written as 

f{x) = n{x)u^{x) (86) 

where n{x) is the density of the corresponding conserved charge in the local rest frame of a fiuid cell. 
( 183|1 and ( 185|) are the equations of motion that need to be solved. An equation of state (EoS) p = p{e) 
is the last equation needed to close the system of equations. It is the only place where the underlying 
dynamics of the system comes into play. In practical applications the only current usually considered 
in heavy ion physics is the net-baryon current j^. Even that current is often negligible at RHIC and 
LHC. 

The equation of state gives us direct information about the QCD phase diagram. This direct link to 
the main goal of the RHIC program is one of the outstanding features of hydrodynamic models. Most 
hydrodynamic computations use an EoS with a first-order phase transition based on the Bag Model. 
However, recent lattice QCD results suggest that the phase transition at low baryon densities is rather 
a smooth crossover |169l 1170] . Since then parameterizations of the EoS from lattice QCD have become 
fashionable. We will describe the recent developments in this area in more detail in Sec. 13.2.31 

3.1.2 Dissipative Corrections 

When we start to include effects of dissipation into relativistic hydrodynamics, we are confronted with 
a rather complicated situation. One of the difficulties is that a naive introduction of viscosities causes 
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the first order theory (i.e. first order in gradients) to exhibit acausahties. Heat might propagate instan- 
taneously because of the parabohc character of the equations. This problem is unique to the relativistic 
generalization of viscous hydrodynamics. In order to avoid this problem, second order terms in heat flow 
and viscosities have to be included in the expression for the entropy |17H I172[ I173[ I174[ I175[ I176[ 1177] , 
but the systematic treatment of these second order terms is difficult. Although there is remarkable 
progress toward the construction of a fully consistent relativistic viscous hydrodynamic theory there 
are still ongoing discussions about the correct formulation of the equations of motion and about the 
appropriate numerical procedures |178] . 

The basic tenet that has to be given up in dissipative hydrodynamics is the assumption of a uniquely 
defined local rest frame. Away from equilibrium the vectors defining the flows of energy, momentum 
and conserved currents might be misaligned. We can still deflne a local rest frame by just choosing a 
velocity u'^(x) in the lab frame. Then the energy-momentum tensor and the conserved current take the 
more general shape 

T^^'ix) = [e{x) + p{x) + U{x)]u^'{x)u''{x) - [p{x) + Il{x)]g^"' + 2W'^^'u''^ + ti^"" , (87) 
f{x) = n{x)u^ + , (88) 

where and are corrections to the flow of conserved charge and energy that are orthogonal to 
and T^^Uy resp., ti^^ (with the orthogonality conditions u^n^'^ = n^^Ui, = 0) is the symmetric, trace-less 
shear stress tensor, and 11 is the bulk pressure. (. . .) around indices indicate symmetrization. Usually 
u'^ is chosen to deflne one of two standard frames: the Eckart frame where the velocity is given by 
the physical flow of net charge (then = 0), or the Landau frame where the velocity is given by the 
energy flow (then = 0). We refer the reader to the article by Muronga and Rischke for further 
details [175] . 

At flrst order the new structures have to be proportional to gradients in the velocity fleld u'^, and 
only three proportionality constants appear in these relations: the shear viscosity t], the bulk viscosity 
( and the heat conductivity k. These transport coefficients are the fundamentally new quantities in 
dissipative hydrodynamics. With the usual definitions the first order relations in the Landau frame are 

m 

n = -cvx (89) 

1' = -^^^'j^ (90) 
e + p T 

vr'^- = 2r/V<''M"> . (91) 

Here g'^ = — (e + p)/nV'^ is the heat fiow in the Landau frame, = {g^"^ — u^u'^)dy is the covariant 
derivative orthogonalized to the fiow vector, T and are temperature and chemical potential for the 
conserved current resp., and (■ ■ ■ ) refers to a symmetrization of indices with the trace subtracted. The 
entropy current receives additional contributions beyond the equilibrium term su^ and one can show 
that with all three transport coefficients positive the entropy is strictly non-decreasing, d^j^S^ > 0. 

At second order many more new parameters, related to relaxation phenomena, appear. Currently, 
most viscous hydrodynamic calculations use the relativistic dissipative equations of motion that were de- 
rived phenomenologically by Israel and Stewart |174] and variations of those, while some use the method 
by Ottinger and Grmela |175[ I176| 1177] , see e.g. |180] . Recently, a second-order viscous hydrodynamics 
from AdS/CFT correspondence was derived |181] . as well as a set of generalized Israel-Stewart equa- 
tions from kinetic theory via Grad's 14-momentum expansion which have several new terms |182] . On 
the other hand however, a stable first-order relativistic dissipative hydrodynamic scheme was proposed 
on the basis of renormalization-group methods |183[ 1184] . The discussion surrounding the appropriate 
(second order) equations of motion is ongoing and a review is beyond the scope of this phenomenological 
overview. We refer the reader to the original articles cited in this subsection for further guidance. 
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In heavy ion physics, most of the focus in viscous hydrodynamics has been on the shear viscosity, 
and its ratio with the entropy density, rj/s. Interesting initial investigations of the effects of bulk 
viscosity have begun |185t 1186] , while heat conductivity still has not been investigated systematically 
in connection with RHIC data. 

3.1.3 Numerical Calculation 

In most ideal hydrodynamic models employed at RHIC numerical computations are carried out in the 
Eulerian formalism. For shock-capturing schemes the SHASTA |187j and RHHLE |188] algorithms 
were developed. The SHASTA algorithm is the most widely spread numerical implementation for 
both ideal and viscous hydrodynamics. On the other hand, the NEXSPHERIO code [189j is based 
on smoothed particle hydrodynamics (SPH) |19Uj . In order to describe the flow of a fluid at high 
energy but low baryon number entropy is taken as the SPH base. Lagrangian hydrodynamics is also 
used in hydro codes for RHIC physics [191J. Lagrangian hydrodynamics has several advantages over 
Eulerian hydrodynamics when applied to ultra-relativistic nuclear collisions. At high energies, the initial 
distribution of the energy is strongly localized in longitudinal direction due to the Lorentz contraction 
of the two nuclei in the lab frame. In order to handle this situation appropriately a fine resolution 
is required in Eulerian hydrodynamics, and computational costs become large. On the other hand, 
in Lagrangian hydrodynamics the grid moves along with the expansion of the fluid. Therefore, one 
can perform the calculation at all stages on those lattice points which were prepared for the initial 
conditions. Another merit of Lagrangian hydrodynamics is the fact that it enables us to derive the 
physical information directly because it traces the flux of the currents. For example, the path of a 
volume element of fluid in the T-^b plane, spanned by temperature and baryon chemical potential, can 
be easily followed. This allows us to discuss directly how the transition between the QGP phase and 
the hadronic phase affects physical phenomena. 

For relativistic viscous hydrodynamics the numerics becomes more complicated in terms of numer- 
ical viscosity and the stability of the calculation. Currently several numerical calculations have been 
implemented and first quantitative comparisons with experimental data become available. Most numer- 
ical algorithms used in viscous hydrodynamics are based on SHASTA |192[ 1193] . but other numerical 
procedures have also been explored, e.g. SPH [194] and the discretization method in Ref. |180] . We 
have only begun to explore the problems related to numerically solving causal relativistic dissipative 
fluid dynamics; see e.g. Ref. |195] for a test procedure applying an algorithm to the Riemann problem. 

3.2 Applications to RHIC Physics 

3.2.1 Hydrodynamics for Heavy Ion Collisions 

Let us first review why RHIC data strongly suggests a long hydrodynamic expansion and cooling phase 
of the bulk matter created in nuclear collisions. We can see the clear signals of bulk collectivity in 
experimental data most convincingly in data on the mean transverse momentum {Pt) as a function 
of observed particle masses and in measurements of elliptic fiow V2 at RHIC. In Fig. |9]we see data on 
{Pt) for several hadron species measured in Au+Au collisions at ^/snn = 130 GeV at RHIC, together 
with a band resulting from the hydrodynamically inspired fit to the tt, K, p, and A data p.96j. The 
observed mass dependence clearly shows that radial fiow exists and that it dominates the motion of 
bulk particles. Small discrepancies between data and theory exist and will be discussed in Sec. HI 

The second remarkable success of hydrodynamics is the explanation of the large elliptic fiow V2 as 
a function of Pt seen at RHIC, and the characteristic dependence on mass, as shown in Fig. |23 To 
realize such strong elliptic fiow in a leading order parton cascade model, unrealistic large cross sections 
are needed p.97] , which clearly supports a hydrodynamic interpretation. The hydrodynamical analyses 
also indicate that thermalization of the matter at RHIC is achieved very rapidly. Most estimates put 
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Figure 9: Mean transverse momentum (Pt) as a function of particle mass in 
central Au+Au collisions at ^/sNN = 130 GeV. Data are from STAR, yellow 
band from a hydrodynamically inspired fit, and the dashed line represents 
the expected behavior from a hydrodynamic fit at the chemical freeze-out 
temperature of 170 MeV at which the radial flow is still too small. Figure 
reprinted from |196j with permission from the American Physical Society. 



the equilibration time between 0.6 and 1.0 fm/c after the collision [3]. f 2 is a wonderful quantity to 
analyze with hydrodynamic models, even deviations from the expected equilibrium values can provide 
information on viscosities |168[ 1198] . 

Hydrodynamics assumes thermal equilibrium and very short mean free paths of particles. We know 
that these assumptions are broken at least at very early times and at the latest times in collisions. They 
might also be broken at the boundaries of the system where densities are smaller, leading to corona 
effects. This means that hydrodynamic calculations have to be combined with some calculations or 
parameterizations of the initial state, and one has to be careful about the treatment of the freeze-out, 
the conversion of hydro fluid cells back into particles. We can further infer from the P^-dependence 
of physical observables such as elliptic flow that thermalization of particles starts to be incomplete at 
high Pt- At RHIC hydrodynamic models generally work very well up to Pt ~ 2 GeV/c. Above that 
threshold novel effects, like quark recombination, come into play at intermediate values of Pt- At the 
largest values, above 6 GeV/c, perturbative QCD production dominates. Of course, the hydrodynamic 
regime comprises about 99% of the particles produced in a typical Au+Au collision at RHIC. 



3.2.2 Initial Conditions 

The hydrodynamic equations of motion need initial conditions for all their dynamic variables, which are 
then evolved forward in time by the solutions. These initial conditions are outside of the framework of 
hydrodynamic models and have to be determined by other means. Physically, they are determined by 
the processes during the initial collision of the nuclei and the approach to equilibrium which is eventually 
reached at a time tq. As we already discussed this equilibration time has been estimated to be rather 
short. However, in practice it should be a parameter since the precise equilibration mechanisms are 
still under debate and calculations of the initial conditions at the start of the equilibrated plasma phase 
have been elusive so far 1200] . 

Historically, parameterized initial conditions for entropy densities (or alternatively the energy den- 
sities) and the net baryon densities have been used |19H I20H I202[ 1203] . In the transverse plane these 
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Figure 10: Initial energy density e(ro, t], x, y) from a Glauber model in longitudinal direction (left panel) 
and in transverse direction (right panel)) for central Au+Au collisions. Figures taken from |191] . 



distributions are usually parameterized based on Glauber-type models of nuclear collisions. In longi- 
tudinal direction often initial distributions inspired by Bjorken's scaling solution are used. Then few 
parameters remain in these initial conditions, such as the maximum values of the energy or entropy 
density, and net baryon density. They are usually fixed by comparison with experimental data on single 
particle rapidity and transverse momentum spectra. 

Let us discuss a straight forward set of energy density and net baryon density distributions that give 
reasonable results. We assume a factorization into a longitudinal profile H{ri) and a transverse profile 
W{x^ y; b) which are given by 

(^{x,y,ri) = e^^^W{x,y;b)H{r]), 
nB{x,y,r]) = nBmi,JV{x,y;b)H{r]) . (92) 

Here the maximum values emax and n^max are parameters, 6 = |b| denotes the impact parameter of the 
collision, and t] is the space-time rapidity. The longitudinal distribution can be parameterized by 

H{r]) = exp [-(Ir^l - Vo)/2al] 9^ - r^o) + ^(r^o - \v\), (93) 

where parameters rjQ and can also be determined by comparison with experimental data on single 
particle distributions. The function W{x, y; b) in the transverse plane is determined by a superposition 
of the density of wounded nucleons — characteristic for soft particle production processes — and the 
density of binary collision — characteristic for hard particle production processes [201] 

W{x,y;b) = w—^ + {l-w)—^-^. (94) 
The density of wounded nucleons is given by 



WN 



T^(b^) (1 - e-^-(^-)'^) + TBihs) (1 - e-^-(^-)'^) (95) 



where = s -|- b, b^ = s — b, and o" ~ 42 mbarn is the total nucleon-nucleon cross section at 
^/sNN = 200 GeV |204j . Ta is the nuclear thickness function of nucleus A, 

Ta{s) = J dzpA{z,s), (96) 
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where the nuclear density pa(^, s) can e.g. be taken to be of Woods-Saxon shape 



PA{r) = po 



1 



(97) 



I _|_ g(r-i?A)/a ■ 



In Eq. fl97|) appropriate parameters a, Ra, and po are 0.54 fm, 6.38 fm and 0.1688 fm ^, respectively 
[204] . On the other hand, the distribution of the number of binary collisions is given by 



The weight factor w can be set to 0.6, consistent with experimental data. Figure [TD] shows examples 
of initial longitudinal and transverse profiles. Obviously the parameters in initial conditions are truly 
additional uncertainties and make quantitative predictions with hydrodynamic models more difficult. 

As a straight forward solution one can choose to set the initial longitudinal flow to Bjorken's scaling 
solution |160] . and one can set the initial transverse flow to zero. This simplest initial flow proflle will 
serve as the basis for all further investigation here. The possibility of an initial transverse flow was 
discussed e.g. by Kolb and Rapp |205] . The initial flow improves the results for P^-spectra and reduces 
the anisotropy. Utilizing a parameterized evolution model it has been pointed out that a Landau-type 
initial condition with complete longitudinal compression and vanishing initial flow is favorable for the 
description of the Hanbury Brown-Twiss (HBT) correlation radii |206[ 1207] . This suggests that HBT 
analyses may be a sensitive tool for the determination of the initial longitudinal flow. As we will discuss 
in detail in the result section, hydrodynamic calculations during the early RHIC years did show very 
bad agreement with experimental data, especially for the ratio of -Rout/-Rside, leading to the notion of 
an HBT puzzle [208] . 

Let us come back to the the apparent early thermalization times found at RHIC. Usually it is said 
that small initial times Tq are needed to describe the elliptic flow data as elliptic flow builds up at 
the earliest stage of expansion when the eccentricity of the flreball is largest [164[ 1165] . However, we 
note that with suitable sets of initial conditions and freeze-out temperatures in fact a larger initial 
proper time is also compatible with data. Luzum and Romatschke show that three very different sets 
of initial and freeze-out temperature (Tj,Tj) — (0.29,0.14) GeV with Tq = 2 fm, (0.36,0.15) GeV with 
To = 1 fm, and (0.43,0.16) GeV with tq = 0.5 fm — provide almost identical differential elliptic flow 
in a viscous hydrodynamic calculation [209] . This suggests that better constraints on initial conditions 
are indispensable to avoid wrong conclusions from comparisons of hydrodynamic calculations with 
experimental data. 

Other approaches for generating realistic initial conditions beyond the Glauber-based parameter- 
izations are available. Color Glass Condensate-inspired initial conditions are becoming increasingly 
popular, see Fig. [TT] They feature larger eccentricities than Glauber-based models which has signifl- 
cant implications on elliptic flow [210] . In that case additional dissipation during the early quark gluon 
plasma stage is needed in order to achieve agreement with experiments [210[ 1211] . Others models are the 
string rope model [212] and the pQCD + saturation model [213] . In the latter the initial time tq is given 
by the inverse of the saturation scale, which is set to a very short r =0.18 (0.10) fm at RHIC (LHC). 
More recently there is a push to implement effects of event-by-event fluctuations in initial conditions. 
In the NEXSPHERIO hydro model these events are created by the event generator NeXus [214] . Fig. 
[T2] shows an example from Ref . [215] . We will discuss the implications of fluctuations in more detail 
later. 

3.2.3 QCD and Hydrodynamics 

Recent lattice QCD calculations show that the phase transition at vanishing or low chemical potential 
is a crossover [169[ I170[ 1216] . However the Bag Model-based equation of state with a flrst-order phase 
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Figure 11: Initial energy density from the Color Glass Condensate in longitudinal direction (left panel) 
and in the transverse plane (right panel) for central Au+Au collisions at ^/snn = 200 GeV. Figures 
reprinted from [210j with permission from Elsevier. 



transition has dominated calculations in the early RHIC years. There are several simulations of lattice 
QCD at high temperature and low density. However, at high net baryon density the investigation with 
lattice QCD becomes difficult. Because of the sign problem, we can not apply the usual Monte Carlo 
methods for finite density lattice QCD. 

A comprehensive study of the dependence of hydrodynamic expansion on the equation of state 
was carried out by Huovinen |2T7] . Figure [13] shows different equations of state used in Ref. [217] . 
Suitable parameters for the initial conditions and freeze-out temperature were chosen for each EoS so 
that the hydrodynamic calculation reproduces the experimental data reasonably. Huovinen found that 
the closest fit to the data of V2 as a function of Pt was obtained for a strong first-order phase transition 
while the results for the lattice-inspired EoS fit the data badly. Similar results were also obtained 
in more detailed analyses |218] . The current discrepancy between lattice results (crossover) and the 
hydrodynamic analysis of data (first order) is unsatisfactory. Dependences on parameters, e.g. fiow in 
the initial conditions, and the freeze-out procedure obscure the observation of the equation of state 
in hydrodynamic calculations. Nevertheless, hydrodynamics offers the most direct tool to investigate 
the equation of state and we have to move to a comprehensive analysis including viscosities, hadronic 
re-interactions and an honest assessment of uncertainties from initial conditions to solve this problem 
in the near future. 

One of the most fascinating topics in connection with the QCD phase diagram is the QCD critical 
point (QCP) and its possible manifestation in relativistic heavy ion collisions |219j . The possible 
existence of this point and its location in the QCD phase diagram have been attracting a great deal 
of interests. Recent studies based on effective theories show a wide range of possible locations of the 
critical point in the phase diagram |220] . In addition, a recent study on the curvature of the critical 
surface in lattice QCD shows that the existence of the QCD critical point is less than certain |221] . A 
solid experimental result seems the best way to settle the question of the existence of the QCP because 
a unanimous conclusion from lattice QCD and effective theories appears to be difficult. For quantitative 
tests of the existence and location of the QCP we need to run hydrodynamics with equations of state 
with a critical point included, and we need to identify appropriate physical observables which show clear 
signatures of a critical point |222[ 1223] . The upcoming lower energy scan at RHIC which will produce 
QCP at higher baryon chemical potential will provide a suitable data sample. 

Relativistic viscous hydrodynamic models need more information besides the equation of state in 
order to run: shear and bulk viscosities, heat conductivity, and relaxation times. On the fiip side one 
should argue that there is a unique chance at RHIC to extract quantitative values for some of these 
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Figure 12: Examples of initial energy density profiles in the transverse plane for central Au+Au collisions 
given by NeXus at mid-rapidity. The energy density is plotted in units of GeV/fm^. Left panel: one 
random event. Right panel: average over 30 random events (similar to the smooth initial conditions 
often used). Figure reprinted from |215j with permission from the American Physical Society. 



transport coefficients from data. Estimates for some transport coefficients are available from lattice 
QCD, finite temperature QCD perturbation theory and effective theories. For strongly coupled quark 
gluon plasma lattice QCD should be a reliable tool. However the determination of transport coefficients 
in lattice QCD with the Green-Kubo formula is not easy. Pioneering results for shear and bulk viscosities 
with large uncertainties have been obtained in |224[ I225j . For weakly coupled QGP perturbative results 
have been obtained by Arnold, Moore and Yaffe [2261 1227[ I228j . Microscopic transport models have 
helped to estimate transport coefficients of a hadron gas |229[l230[l231j and a parton gas [232[l233j . The 
results seem to indicate that t]/ s from a hadron gas is larger than the KSS bound l/(47r) |181j which is 
proposed by the AdS/CFT correspondence. On the other hand, shear viscosities from radiative parton 
transport can be surprisingly lower than leading order perturbative results and close to the AdS/CFT 
bound. We have accumulated a few predictions in Fig. [3H Even though there is still little trust in 
quantitative numbers the preliminary conclusion is that the origin of the perfect fluid signatures at 
RHIC should be in the partonic phase |231j . 

There are also attempts to compare hydrodynamic calculations to partonic transport models. Early 
parton cascades with leading order dynamics |234] had suggested that large cross sections are needed 
to account for the elliptic flow observed at RHIC. Thermalization of the parton matter in the initial 
stage would be extremely slow. More recent parton cascades including radiative corrections have led 
to results that are compatible with large elliptic flow and rapid thermalization [2351 12361 1237] . More 
studies comparing hydrodynamics and transport models in detail are needed as they provide necessary 
mutual checks. 

3.2.4 The Freeze-out Process 

Currently two separate freeze-out processes are believed to occur in heavy ion collisions. One is the 
chemical freeze-out at which the ratios of hadrons are flxed, and the other is the kinetic freeze-out at 
which the particles stop to interact. Chemical freeze-out temperature and chemical potentials are ex- 
tracted with the statistical model on the basis of the grand-canonical formalism. Surprisingly, statistical 
models are in excellent agreement with experimental data for hadron ratios in a wide range of collision 
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Figure 13: (a) The entropy density s divided by T^, and (b) the energy density e divided by T"^ both 
as functions of temperature; (c) the pressure p, and (d) the velocity of sound squared cj. = dp/de both 
as functions of e. In all four panels four equations of state are shown: (i) qp = quasiparticle EoS, (ii) 
Q = ideal parton gas with first order phase transition, (iii) H = hadron resonance gas, and (iv) T = 
ansatz with crossover. Figure reprinted from |217] with permission from Elsevier. 

energy from AGS to RHIC |238j . From the values of chemical freeze-out temperatures around 170 MeV 
we can infer that the chemical freeze-out process takes place just after the QCD phase transition [238]. 
At the kinetic freeze-out temperature the mean-free path of particles has grown to be of the order of 
the system size and a hydrodynamic description which requires zero mean-free path is clearly no longer 
applicable. A first naive guess for the kinetic freeze-out temperature Tf is the pion mass (~ 140 MeV). 
Practically the value of the kinetic freeze-out temperature in a hydrodynamic model is determined from 
comparison with data for P^-spectra. 

The task at the end of a hydro calculation is the population of fluid cells of a given temperature 
and flow with particles. For calculations of single particle spectra, the simple assumption of a sudden 
freeze-out process at a certain proper time for each fluid cell is sufficient, neglecting the reverse process 
from particles to the hydrodynamic medium. Under this assumption the Cooper-Frye formula |239] is 
widely used. Here we start from a simple case |240] : Suppose that a number of particles N{t) exists in 
the enclosed volume Q that is bounded by a closed surface S{t) at time. N{t) is given by 



where n(r, r) is particle number density. At time r + 6t, the number of particle N{t) has changed to 



where the volume has changed to + dQ. Utilizing current conservation, ^ + V ■ j (j is the current 
of particles), Eq. f llOOp is rewritten as 




(99) 




(100) 




(101) 
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Table 2: Calculations with relativistic ideal hydrodynamical models at RHIC 
energies. See text for acronyms. 



Ref. 




IC 


EoS 


Hydro 
Exp. 


hadronization & FO 


observables 


Kolb 


164J 


Glauber 


1st order 


2d 


KF 


elliptic flow 


Huovinen \W5\ 


Glauber 


1st order 


2d 


KF 


radial, elliptic flow 


Kolb |201j 


Glauber 


1st order 


2d 


KF 


centrality dependence 














of multiplicity and 














flow 


Hirano |2()2] 


Glauber 


1st order 


3d 


CF and KF 


flow and HBT 


Teaney [203J 


Glauber 


1st order 


2d 


hadron cascade 


flow, QGP signature 


Nonaka [19 Ij 


Glauber 


1st order 


3d 


hadron cascade 


SDectra flow 


Kolb |205j 


Glauber, 


1st order 


2d 


CF and KF 


Pt spectra, V2 






initial Vt 










Eskola |213j 


pQCD + 


1st order 


cylindrical 


KF 


low and high Pt spec- 






saturation 




symmetry 




tra 


Hirano |210j 


CGC 


1st order 


3d, hy- 


CF and KF 












dro+jet 






Hirano [211] 


CGC 


1st order 


3d 


hadron cascade 




Andrade [m] 


NEXUS 


1st order 


3d 


KF 


f 9 


Hirano |287j 


CGC + 


1st order 


3d 


hadron cascade 


V2 






fluctuation 










Huovinen [217] 


Glauber 


comparison 


2d 


KF 


hadron spectra, f 2 


Socolowski j215i| 


NeXus 


1st order 


3d 


CE 


HBT 


Hirano |267j 


Glauber 


1st order 


3d, hy- 


CF and KF 


hadron spectra, V2 










dro+jet 







where n is the normal vector of surface element d'^s and d( is the distance between the surface of ^{t) 
and that of ^{t) + dQ. In Eq. fllOip . dN/dr is the number of particles which cross the surface S{t) 
during dr. Then the total number of particles through the hypersurface S, which is the set of surfaces 
{Sir)}, is 

N = [ fda^, (102) 
where j° = n, dao = d^r, da = drd'^s n. If we write 



E (27r)3exp[(P,M--/i;)/T;]±l ^ ^ 

for the current in Eq. f ll02p . we obtain the Cooper- Frye formula |239] 

d^P V (2^)' exp[(P,n^ - ^,f)lTf] ± 1' ^''^^ 

where qh is a degeneracy factor of hadrons and Tf and /i/ are the freeze-out temperature and chemical 
potential. In other words we obtain da^j, by estimating the normal vector on the freeze-out hypersurface 
E. Using Eq. f ll04p . we can then calculate all particle distributions after freeze-out. 

More realistic models have been investigated. One of them is the Continuous Emission Model (GEM) 
in which particles are emitted continuously from the whole expanding volume of the system at different 
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Table 3: Calculations with relativistic viscous hydrodynamical models at RHIC 
energies. See text for acronyms. 



Ref. 


IC 


EoS 


Hydro 
Exp. 


hadronization & 
FO 


observables 


Song[192j 


Glauber 


EoS de- 
pendence 


2d, I-S 


KF (direct vr) 




DuslingflSD] 


Glauber 


ideal QGP 

gas 


2d, 0-G 


KF, viscous cor- 
rections 




RomatschkepSS] 


Glauber 


semi- 
realistic 
EoS 12471 


2d, I-S 


KF 


multiplicity, V2 


Luzum|2U9j 


Glauber 
and CGC 


semi- 
realistic 
EoS 1247] 


2d, con- 
formal 
relativistic 
viscous 
hydro [181] 


KF, viscous cor- 
rections 


multiplicity, vt, f2 


(Jhaudhuri|193j 


Glauber 


lattice 

+HRG 

EoS 


2d, I-S 


KF 


V2 



temperatures and different times |215] . In the early days of hydrodynamics only kinetic freeze-out 
was implemented. Indeed, at lower collision energies such as at AGS and SPS, the differences between 
chemical freeze-out and kinetic freeze-out points are not large. However, at RHIC a significant difference 
between kinetic freeze-out temperatures from hydro-inspired models and the chemical freeze-out from 
the statistical model appears |241] . This phenomenon also manifests itself through the failure to get 
the correct absolute normalization of some P^-spectra, e.g. the proton in hydrodynamic calculations, 
with only a kinetic freeze-out |242j . Hence a consistent modeling of separate freeze-out processes via 
modified equations of state was introduced [202i I243i I244i 1245] . 

It turns out that some experimental data is still not understood in a satisfying way even with two 
separate freeze-out procedures. For example, mean transverse momentum (Pt) as a function of particle 
mass does deviate from the linear scaling law, which suggests significant final state interactions in the 
hadronic phase |191] . To explain these effect, and to account for the apparently large viscosities in 
the hadronic phase, as discussed before, hydro+cascade hybrid models were introduced. They use a 
hydrodynamic computation of the expansion and cooling of hot QCD bulk matter, and then couple the 
output consistently to a hadron-based transport model for the final-state interactions. Pioneering work 
on hydro+cascade hybrid models was done by Bass et al. )i246j using UrQMD. Similar investigations 
were carried out in Refs. [203^ 1211] . The improvements introduced by these hybrid models are discussed 
in more detail in the results section. 

3.3 New Developments 

Through its success at RHIC hydrodynamics has positioned itself as one of the most useful phenomeno- 
logical models for heavy ion collisions. A large amount of studies working with hydrodynamics have 
been carried out in the RHIC era. Historically, most of them have been using ideal hydrodynamics, but 
obviously viscous hydrodynamics will be a major focus point for the near future, as both mathematical 
and numerical issues are settled. The thorough vetting of hydrodynamic models has also become a 
top issue. Different existing codes have to be tested against each other using identical parameters and 
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Figure 14: Centrality dependence of pseudora- 
pidity distributions of charged particles from 
the hydro+UrQMD hybrid model [191 j com- 
pared wit PHOBOS data for ^IJ^ = 200 GeV 
Au+Au collisions |270j . Impact parameters in 
the calculation with the corresponding central- 
ity bins in parentheses are 2.4 (0-6%), 4.5 (6- 
10%), 6.3 (10-15%), and 7.9 fm (25-35%). Fig- 
ure taken from [191] . 
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Figure 15: P^-spectra for 7r+, and p at cen- 
tral y/SNN = 200 GeV Au+Au collisions with 
PHENIX data |271j . The points are not renor- 
malized. Figure taken from |191] . 



initial conditions. One such effort is under way in the "Theory-Experiment Collaboration for Hot QCD 
Matter" (TECHQM) [2l8] . 

Hydrodynamic models can be easily adopted for descriptions of the entire process of heavy ion 
collisions that include initial collisions, thermalization, bulk dynamics, hadronization, freeze-out, final 
state interactions and even hard and electromagnetic probes. Such comprehensive approaches are called 
multi- module modeling. There are ongoing efforts to construction such models [191^ 1211^ 1249^ 1250] . The 
coupling of hydrodynamics and hadronic transport together to form a hybrid model is only the first 
step. Event generators for the initial state will soon become standard. It has also become necessary to 
conduct jet quenching studies using realistic input from hydrodynamically evolved fireballs |152j . 

We have compiled a list of relativistic ideal hydrodynamic models which are mentioned in this work 
in Tab. |2l We reference each work and quote the initial conditions, freeze-out treatment, and equation 
of state that have been used, and the observables that have been calculated. The acronyms KF, CF 
and CE stand for kinetic freeze-out, chemical freeze-out and continuous emission, respectively. 

For relativistic viscous hydrodynamics, the number of phenomenological studies is much smaller. 
Quantitative discussions comparing to RHIC data have been carried out, but need to be weighed with 
caution. Many points need further study. In Tab.|3l we list the relativistic viscous hydrodynamic studies 
which compare results to experimental data from RHIC. Here I-S (G-0) stands for relativistic viscous 
hydrodynamics proposed by Israel and Stewart |172l 1173^ 1174] (Grmela and Ottinger [175^ 1176^ I177j ). 
The major difference between ideal hydrodynamics (Tab. [2]) and viscous hydrodynamics (Tab. [3]) is that 
the former all agree on the set of equations to be solved, while the latter solve a variety of different second 
order schemes. While all of those should only exhibit small deviations from first order hydrodynamics, 
this strengthens the point that we have to apply due caution when analyzing RHIC data. 
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Figure 16: Pr-spectra of 7r+ from |191j . Shown 
are the spectra for hydrodynamics at the 
switching temperature of 150 MeV, result from 
hydro + decay only, and hydro + full UrQMD, 
in central collision. Figure taken from |191] . 
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Figure 17: Mean Pt as a function of mass. Fig- 
ure taken from [19 Ij. 



3.4 Hadronization and Quark Recombination 

In hydrodynamic models the hadronization process from the QGP phase to the hadron phase is naturally 
encoded in the equation of state. As we have mentioned above hydrodynamic models generally do well 
for particles below Pt ^ 2 GeV/c at RHIC. Between that point and the perturbative region above ~ 6 
GeV/c an intermediate Pt region has been found. It exhibits features from both domains, without 
fitting exclusively in any of the two categories. E.g. we find large elliptic flow and baryon over meson 
ratios that are reminiscent of hydrodynamics and soft physics, and incompatible with jet quenching and 
fragmentation. On the other hand the elliptic flow is not large enough to follow the ideal hydrodynamic 
predictions, and we find dihadron correlations at intermediate Pt that exhibit clear jet-like structures. 
Simple interpolations (hydro+jet models) do not explain all features, e.g. the much smaller elliptic flow 
of mesons compared to protons. In order to explain hadron production at RHIC, and in order to 
understand the breakdown of the pQCD or hydrodynamic approach it is crucial to understand this 
kinematic region at intermediate Pt- 

Quark recombination or coalescence is the best candidate to explain a large amount of experimental 
data at intermediate Pt- Recombination models assume a universal phase space distributions of quarks 
at hadronization. Quarks turn into baryons, qqq — )■ B, and mesons, qq ^ M described either by using 
instantaneous projections of quark states onto hadron states |25H \'25'2\ \25'S\ I254[ 1255] . or a dynamical 
coalescence process with finite width hadrons governed by rate equations |256j . Note that usually only 
the valence quarks of the hadron are taken into account although generalizations have been worked out 
[255] . 

The original instantaneous projection models explicitly preserve only three components of the 
energy-momentum four-vector in the underlying 2 — )■ 1 and 3 — )■ 1 processes. The yield of mesons 
can be expressed through the convolution of the Wigner function Wab for parton pairs a, b and the 
Wigner function encoding the meson wave function 

dNM f d^R f d^qdh- / r P r P \ , , , , 

7^"'- (R + 2. 2 + q; R - 2 , 2 - qj t«(r, q) . (105) 
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Figure 18: Spectra from hydrodynamics (dot- 
ted), pQCD (dash-dotted), and their sum (sohd 
hne) for vr", K~ , and p in Au+Au colhsions at 
^snn = 200 GeV for impact parameter b = 2.0 
fm compared to data from PHENIX [272] . Fig- 
ure reprinted from |267j with permission from 
the American Physical Society. 
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Figure 19: Ratios Np/N^^- and Nk-IN^,^ as a 
function of Pt in Au + Au colhsions as in Fig. 
la compared to data from PHENIX [272]. Fig- 
ure reprinted from |267j with permission from 
the American Physical Society. 



The quark Wigner functions are usually approximated by classical phase space distributions. Hadron 
spectra at intermediate Pt are described well by considering a factorization into thermal quark distri- 
butions [252] 

Wah (ri, pi; Y2, P2) = fa iji, Pi) fb (1-2, P2) • (106) 

Correlations between quarks can be introduced to model correlations found between hadrons |257] 
without interfering with the excellent description of spectra and hadron ratios. 

Dynamical models, like the resonance recombination model |256[ I258[ 1259] solve systems of rate 
or Boltzmann equations for the underlying 2 — )■ 1 and 3 — ?■ 1 processes. In the case of resonance 
recombination mesons and baryons are technically treated as resonances with widths Pj^, and the inverse 
processes (the "decays" of mesons and baryons into quarks) are taken into account to achieve detailed 
balance. Resonance recombination has the advantage that is can deliver an equilibrated hadron phase 
from an equilibrated quark phase due to energy conservation and detailed balance. This has opened 
the possibility to reconcile quark recombination with hydrodynamics and equilibrium physics at low 
momenta [259]. 

At large momenta contact can be made with jet fragmentation by rewriting fragmentation functions 
as a process of valence quark coalescence in a suitably defined jet shower [260[ I26H 1262] . If S{p) is the 
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Figure 20: Top panel: inclusive Py-spectrum of charged hadrons in central Au + Au collisions at 
y/sNN = 200 GeV. Contributions from a recombination model and perturbative liadron production 
are shown compared to experimental data from PHENIX j273j . Bottom panel: ratio of protons 
to positive pions as a function of Pt for the same Au+Au collisions. The region below 4 GeV/c 
is dominated by recombination, the region above 6 GeV/c by fragmentation of hadrons from jets. 
Figure taken from |251] . 



distribution of quarks in a jet before hadronization and T{p) is the distribution of the underlying event 
(partially thermalized in heavy ion collisions) recombination will be applied to the total parton phase 
space S{p) + T{p). For mesons this will lead to the following possible combinations: SS which resembles 
the fragmentation process within a jet, TS which is a novel "soft-hard" hadronization process, and TT 
which corresponds to the usual picture of quark recombination. Shower distributions S{p) can be fitted 
such that SS and SSS recombination fit the known fragmentation functions for the respective mesons 
and baryons. 

The strength of the quark recombination picture is the predictive power coming from explaining all 
measured hadron spectra at intermediate Pt with one parameterization of the quark phase at hadroniza- 
tion. It has been shown that at low momenta resonance recombination is compatible with hydrodynam- 
ics and kinetic equilibrium |259] . but on the other hand, because a thermalized state does not retain 
memories of a previous time evolution, all phenomena in the equilibrated region should be explainable 
by hydrodynamics. This includes the quark number and kinetic energy scaling observed at RHIC at low 
momenta |259j . The possibility of including quark recombination explicitly into hydrodynamic model 
has been studied in [263] . 

The quark number scaling law is a signature feature of quark recombination. For a quark phase with 
elliptic flow V2{Pt) at the time of hadronization simple instantaneous recombination models predict 



where n is the number of valence quarks. This scaling law describes a key feature of experimental 
data at intermediate Pt rather accurately. We make two comments here: (i) the general shape of 
V2 at intermediate momenta suggests that in the kinematic range under consideration {Pt between 
1.5 and 5 GeV/c) hadrons are no longer equilibrated or at least there are large viscous corrections to 
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Figure 21: Raa for neutral pions in Au+Au and Cu+Cu collisions at top RHIC energy 
and for different centralities, calculated in the GLV model. Finalized PHENIX data is now 
available in |275] . Figure reprinted from |277] with permission from Elsevier. 



hydrodynamics; (ii) at lower momentum the scaling seems to work rather with kinetic energy instead 
of Pt- This is a rather accidental feature of V2 for equilibrium hydrodynamics and, as already brought 
up above, can not directly be attributed to quark recombination |259j . 

We discuss the situation from the experimental point of view further in the result section |H More 
profound detours into quark recombination models are beyond the scope of this review and we refer the 
reader to several review articles on quark recombination for further study [2641 12651 1266] . 



4 Interpretation of Experimental Data from RHIC 

We now proceed to discuss some key experimental results from the first decade of running of the 
Relativistic Heavy Ion Collider. Many of those results can be understood, at least qualitatively, within 
the framework of perturbative QCD and hydrodynamics. We will also review some attempts to extract 
quantitative statements about fireball parameters and properties of quark gluon plasma. We proceed 
from the simplest observables, single particle yields and spectra to more intricate ones and try to take 
a comprehensive view that includes both bulk and hard probe particles. 
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Figure 22: Raa for several centrality bins for Au+Au collisions at top RHIC energy. Cal- 
culations using ASW energy loss with reweighted (dashed) and non-reweighted quenching 
weights (solid lines) are compared to data from PHENIX for charged hadrons (closed squares) 
and neutral pions (open squares) [H I278j . and from STAR for charge hadrons P]. Figure 
reprinted from |123j with kind permission from Springer Science + Business Media. 



4.1 Particle Multiplicities and Single Particle Spectra 

Hydrodynamic models by default deliver the main qualitative features of bulk spectra measured in 
heavy ion collisions at RHIC: hadrons at low transverse momentum Pt are thermalized and exhibit 
a relativistic outward collective flow. Nevertheless it is a challenge to quantitatively describe details, 
like the subtle differences seen between hadron species. If we want to describe bulk observables with 
hydrodynamics or hybrid models based on hydrodynamics the first goal is the determination of realistic 
initial conditions. Even if some theoretical modeling of the initial state is available, there are usually 
a few parameters that are fitted to single particle P^-spectra and rapidity distributions in central 
collisions. 

In Figs. [TH and [15] we show single particle spectra from a hybrid model of (3+l)-dimensional ideal, 
relativistic hydrodynamics and UrQMD |191] . In this model final state interactions at the end of the 
fireball life time are taken into account by connecting the hydrodynamic phase to the hadron-based 
event generator. For the determination of initial conditions, this hybrid model works the same way 
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Figure 23: Left panel: Raa of neutral pions for central Au+Au collisions at top RHIC ener- 
gies using LO and NLO accuracy hard processes and medium-modified fragmentation func- 
tions from the Higher Twist formalism compared to data from PHENIX. Finalized PHENIX 
data available in |275j . Right panel: x^-fit of eo to RHIC data on Raa and Iaa yielding a 
consistent value for both observables. Figures reprinted from |279] with permission from the 
American Physical Society. 



as pure hydrodynamics. Practically, the parameters to be determined are the maximum values of 
energy density and baryon number density which are fixed from comparison with experimental data of 
pseudorapidity distributions and Py-spectra in central collisions. The centrality dependence of spectra 
is then settled and reproduces the data well as shown in Fig. [HI 

Figure [13 compares the P^-spectra for vr"'", K and p in central Au + Au collisions. Because the 
hybrid model correctly treats chemical and kinetic freeze-out processes, taking into account the different 
cross sections of hadrons, the absolute value of the proton spectra shows good agreement with the 
experimental data. This is one of the improvements that hybrid models with hadronic transport offer 
over pure hydrodynamic model with only one (thermal) freeze-out process. Such pure hydro calculations 
usually fail to explain proton spectra. Clearly, taking into account the correct freeze-out for each particle 
species separately is important and hybrid models with hadronic transport deliver that. 

Let us now investigate the detailed effect of resonance decays and hadronic rescattering on the shape 
of the momentum spectra. Figure [16] compares the Pp-spectrum for vr"*" from hydrodynamics only, i.e. 
at the switching temperature T^^ = 150 MeV of the hybrid model (solid line), to the spectrum after 
resonance decays have been taken into account in addition (open symbols). We also show the result if 
the full UrQMD transport is run on the hydro result (solid circles). The difference between the solid line 
and open symbols directly quantifies the effect of resonance decays on the spectrum. They obviously 
increase the yield of pions, most dominantly in the low momentum region P^ < 1 GeV/c. Furthermore, 
the difference between open and solid symbols quantifies the effect of hadronic rescattering. Pions with 
Pt > 1 GeV lose momentum resulting in a steeper slope. In other words hadronic re-interactions cool 
the spectra. However, they do it selectively with the appropriate cooling rates for different hadron 
species. 

To strengthen this point we show the mean transverse momentum (Pt) as a function of hadron 
mass from |191j in Fig. [T71 As before we compare results at the switching temperature Ts„ = 150 MeV, 
corrected for hadronic decays (open symbols) to the results from hydro plus full UrQMD (solid symbols). 
In the former case (Pr) follows a straight line, as expected from a hydrodynamic expansion. It is a 
natural consequence of collective flow. If hadronic rescattering is taken into account (Pr) does no longer 



50 



1 

0.8 



<0.6 

< 



0.4 
0.2 

0.8 

<0.6 

< 



0.4 
0.2 






• PHENIX 0-5% 
- AMY, b = 2.4 fm, = 0.33 

~ ~ ri 1 , D = Z.4 rm, (J = 1 .joe v /nn, c^^^ = u.z 
. _ . ASW, b = 2.4 fm, K = 3.6 


' 1 ' 1 

f J . 




T T T T I 


— 1 — 1 — 1 — 1 — 1 — 1 — 1— — 1 — 1 — 1 


1 — 1 — 




■ ■ PHENIX 20 - 30% 

- AMY, b = 7.5 fm, a = 0.33 

- — HT, b = 7.5 fm, ^ J, = 1 .5GeV7fm, c^^^ = 0.2 
. ■ — ■ ASW, b = 7.5 fm. K = 3.6 










: \ 

1 . 1 . 1 . 1 


r I ] 

1 


1 



10 



12 14 
P.J. (GeV/c) 



16 



18 



20 





1.25 




1.2 




1.15 




1.1 












1.05 






< 


1 




0.95 




0.9 




0.85 



■ ■ AMY, P.J, = 10 GeV/c 

AMY, P.J, = 15 GeV/c 

• .HT, p.j.= lOGeV/c 

HT,p,j.= 15 GeV/c 

A A ASW, P.J. = 10 GeV/c 




0.5 



(|) (rad) 



1.5 



Figure 24: Left panels: Raa as function of Pt for central (top) and mid-central (bottom) 
collisions calculated from the ASW, Higher Twist and AMY energy loss models. The single 
parameter in each model has been fitted to describe the data by PHENIX (finalized data 
available in |275j ). Right panel: Raa as function of azimuthal angle normalized by the 
0- integrated value for two different values of Pt- Again all three energy loss models are 
shown. Figures taken from |152] . 



follow the straight line. The average momentum of pions is actually reduced by hadronic rescattering 
(they act as a heat bath in the collective expansion), whereas protons pick up additional transverse 
momentum in the hadronic phase. Data by the STAR Collaboration is shown as well (solid triangles). 
The proper treatment of hadronic final state interactions significantly improves the agreement of the 
calculation with the data. 

Generally, P^-spectra from hydrodynamic calculations show good agreement with experimental data 
up to Pt ~ 2 GeV/c. For larger Pp discrepancies between hydrodynamic results and experimental data 
appear, which suggest the approximate limits of applicability for hydrodynamics. We can see these 
discrepancies between hydrodynamic results and experimental data in spectra in Fig. [18] and in elliptic 
flow in Fig. |27] (which we discuss in more detail later). The general features of extended spectra can 
usually be described in hydro plus jet models, as e.g. shown in Ref. [267]. 

Figure [18] shows Pp-spectra for vr", K~ and p in Au + Au collisions at ^/snn = 200 GeV. The dotted 
lines indicate ideal hydrodynamics. These spectra exhibit a Boltzmann-like shape (~ exp{—E/T)) 
which is locally boosted by a radial flow velocity (m^) through E — )• p^u'^. This is the gross feature 
of all P^-spectra in hydrodynamic models and comparison with data suggests that thermalization is 
reached. The figure also shows the contributions from hard (pQCD) hadron production (dash-dotted 
lines), whose spectra exhibit power-law behavior. Because of the exponential suppression of the bulk 
perturbative hadron production dominates at large enough Pp. The sum of both hydro and perturbative 
("jet") contributions describes the data quite well. The value of Pt where the transition from the soft 
to the hard component takes place depends on particle species. In the case of pions it happens between 
Pt ~ 1 and P^ ~ 2 GeV/c in this particular calculation. For protons the transition happens gradually 
in the range 2 < Py < 5 GeV/c, but certainly at higher momentum than for pions. The crossing point 
is determined by two facts: one is radial flow which pushes the soft components toward high Pp. The 
other is the jet quenching mechanism which suppresses contributions from jets. However, the simple 
hydro-|-jet model does not explain some striking features of the systematics of different particle species. 

The differences in the transition points for different particles come from the different (mass depen- 
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Figure 25: V2{Pt) for pions, kaons and protons 
produced in minimum-bias collisions at RHIC 
from an ideal hydro calculation by Huovinen 
et al. |165j compared to data from PHENIX 
[280] ■ Figure reprinted from |280j with permis- 
sion from Elsevier. 
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Figure 26: V2 as a function of pseudorapid- 
ity 7] for charged hadrons in Au+Au collisions 
at -sJsnn = 130 GeV. The solid and dashed 
lines correspond to the hydrodynamic calcula- 
tions with partial chemical equilibrium (PCE) 
and full chemical equilibrium (CE) respectively. 
Data is taken from STAR [M] and PHOBOS 
Figure reprinted from |202] with permis- 
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dent) effects of radial flow, and from the different yields of these particles in the fragmentation process 
(which disfavors heavier mesons and baryons). The effects of different crossing points are most clearly 
seen in hadron ratios as a function of Pt- In Fig. [19] we show the hadron ratios p/vr" and K~ /tt~ as a 
function of Pt from the hydro-|- jet model advocated in \2E7\ . A crucial point in this figure is the fact 
that the ratio of p/vr" is as large as unity at ~ 3 GeV/c, which can not be understood from pQCD. 
This anomalously large yield of protons at intermediate Pt was the first indication of the so-called 
baryon puzzle. 

It was the baryon puzzle that gave birth to quark recombination models in the RHIC era. Most 
of the experimental evidence for recombination comes from elliptic flow measurements, but the large 
baryon/meson ratios were essential to highlight the necessity for a mechanism that is able to push the 
region of soft physics farther out for baryons than for mesons. One crucial event was the advent of first 
data on the meson. In a hydro-|-jet model 0s essentially behave like protons, because of their mass, 
while in data they exhibit the universal behavior of other mesons (pions, kaons) at intermediate Pt 
[252112681 1269]. 

Figure [20] shows the Pr-spectrum of charged particles and the p/vr^ ratio as a function of Pt from 
the model advocated in [251] which includes quark recombination and hadrons from pQCD ( "fragmen- 
tation"). The proper hydrodynamic region below 2 GeV/c has been omitted. From comparison with 
experimental data from STAR one can conclude that the region below 4 GeV/c is actually dominated 
by quark recombination while the pure pQCD region free of soft or bulk hadron production only starts 
beyond 6 GeV/c. Similar conclusions have been reached in other recombination + jet models. We will 
find even stronger evidence when we discuss elliptic flow in the next subsection, where recombination 
shows clearly visible and experimentally tested differences compared with both hydrodynamics and 
perturbative hadron production. 

Beyond 6 GeV/c in transverse momentum hadron spectra are dominated by perturbative production 
that includes energy loss through final state interactions with the medium and fragmentation. Since 
the interesting observable is the suppression with respect to high-P^ jets and hadrons in the vacuum 
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Figure 27: Elliptic flow V2{Pt) for vr^, K^, 
and p in Au+Au collisions at ^/s^n = 200 
GeV at impact parameter h = 5.5 fm from the 
hydro+jet model |267] . Figure reprinted from 
|267j with permission from the American Phys- 
ical Society 
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Figure 28: Upper panel: V2 for p and vr^ from 
the recombination plus fragmentation model in 
(251] , Lower panel: V2 for and A + A from 
the same model. Data shown is from PHENIX 
[282j and STAR [283j. Figure taken from j25T] . 



the best way to analyze the data is by considering the ratio 

of yields in A+A vs p + p collisions. This nuclear modification factor is similar to what we had defined 
in Eq. ( l38l) for d+A collisions. (A'^cou) is the average number of binary nucleon-nucleon collisions that 
we expect for a given centrality bin. It is usually determined from Glauber-type model calculations. An 
incoherent superposition of nucleon-nucleon collisions would lead to Raa = 1 modulo isospin effects. 
We have already discussed in Sec. 12.2.41 how initial state nuclear effects lead to modest deviations from 
unity. 

As had been predicted, first RHIC data at high Pt revealed a huge suppression of hadrons [H [21 
I274[ I275[ 1276] . This strong jet quenching is one of the pillars on which a qualitative argument for the 
creation of quark gluon plasma at RHIC energies rests. Fig. [21] shows calculations by Vitev within the 
GLV model for the Raa of neutral pions in Au+Au and Cu+Cu at RHIC energies of ^/sNN = 200 
GeV for different centralities [277] . The calculation also takes into account nuclear effects in the initial 
state. The quenching is as large as a factor 5 in central Au+Au collisions as indicated by the data from 
PHENIX. Vitev connects the quenching strength g to a local gluon density dNg/dy in the medium and 
he finds dNg/dy ~ 800 — 1175 for central Au+Au collisions. 

Figure [22] shows a systematic study of the centrality dependence of the Raa of charged hadrons 
in Au+Au collisions carried out by Dainese, Loizides and Paic using an implementation of the ASW 
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Figure 29: PT-integrated elliptic flow for charged hadrons at midrapidity {\ri\ < 1) from Au + Au 
collisions at RHIC top energy as a function of the number of participating nucleons iVpart- Thin 
lines show the predictions from ideal fluid dynamics only with a freeze-out temperature Tjec = 100 
MeV while thick lines indicate the results from the hydro+JAM hybrid model. For both models two 
different initial conditions, CGC and BGK, are used. Figure reprinted from |211] with permission 
from Elsevier. 



energy loss model |123] . They model the centrality dependence by scaling not only the production 
of jets but also the local q by the density of binary nucleon-nucleon collisions, g(x) = A;ncoii. This is 
different from other approaches that usually scale soft particle densities in the transverse plane by the 
density of nucleon participants npart- The centrality dependence of the data from PHENIX and STAR 
is reproduced well by just one adjustable parameter, the normalization k. So far, no conclusive picture 
has emerged which scaling of q is favored by data. The fit to data is aided by the large theoretical 
uncertainty that comes from the difference between the reweighted and non-reweighted versions of the 
ASW quenching weights that have been discussed earlier. Those two results define the grey bands in Fig. 
[22l The average quenching found by the authors within the ASW formalism is g ~ 15 GeV^/fm. This 
value, like others obtained with the ASW model, are rather large compared to those from alternative 
energy loss calculations as we will see below. 

Zhang et al. have presented one of the first studies of jet quenching using next-to-leading order hard 
processes [279] . The left panel of Fig. [23] compares their results for the Raa of 7r° in central Au-|-Au 
collisions at LO and NLO accuracy with data from PHENIX. They use medium-modified fragmentation 
functions inspired by the Higher Twist formalism that are rescaled by the average energy loss of a parton, 
and then use these instead of NLO vacuum fragmentation functions in a next-to-leading order code for 
hadron production. They parameterize the energy loss through an integral over a transverse profile of 
an expanding fireball along the path of the parton times a normalization parameter eo that characterizes 
the stopping power parameter of the medium. The left panel of Fig. [23] shows their result for three 
different values of e that define an approximate uncertainty band. Both the LO and NLO calculations 
can describe the data if eo is a fit parameter. If eo is kept fixed the NLO calculation shows stronger 
quenching due to the larger ratio of gluon to quark jets which couple more strongly to the medium. 
The right panel of Fig. [23] shows a x^-fit of to Raa between 4 and 20 GeV/c and to Iaa (discussed 
below), yielding consistent values in a range eo = 1.6 . . . 2.1 GeV/c |279j . 

Fig. [2^ shows results from a comparative study by Bass et al. Jets are propagated through a medium 
described by hydrodynamics, using three different schemes for energy loss: ASW, Higher Twist, and 



54 



0.1 
0.08 
0.06 
0.04 
0.02 



0.1 
0.08 
0.06 
0.04 
0.02 



(b) 15-25% b=6.3fm 



0.12 
0.1 
c 0.08 
> 0.06 
0.04 
0.02 



(a) 3-1 5% b=4.0f m • hydro+cascade 
^ ' "T =100MeV 



Tf'=169MeV 

dec 

PHOBOS 



oi 



1 • 



i • 



•-• • • • • 



(c) 25-50% b=8.5fm 



-5 -4 -3 -2 -1 1 2 3 4 5 



0.1 
0.08 
0.06 
0.04 
0.02 



0.1 
0.08 
0.06 
* 0.04 
0.02 

ot 

0.12 
0.1 
c 0.08 
> 0.06 
0.04 
0.02 



(a) 3-15% b=4.0fm • hydro+cascade 
^ ' t.uiiii ^ =ioOMeV 



"=1 69MeV 

dec 

PHOBOS 



■ 



(b) 15-25% b=6.3fm 





-5 -4 -3 -2 -1 1 2 3 4 5 



Figure 30: Left panel: The pseudorapidity dependence of V2 for charged hadrons in (a) central (3- 
15%), (b) semi-central (15-25%), and (c) peripheral (25-50%) Au + Au collisions at y/sj^ = 200 GeV. 
The corresponding impact parameters are b = 4.0, 6.3, and 8.5 fm resp. The hydrodynamic evolution 
is initialized with modified BGK initial conditions. The lines show the predictions from ideal fluid 
dynamics only for freeze-out temperatures of T^ec = 100 MeV (solid blue) and Tjec = 169 MeV (dashed 
green). Red circles show the corresponding results from the hydro+JAM hybrid model. Right panel: 
Same as left panel, except for using CGC instead of BGK initial conditions. Figures reprinted from 
[211] with permission from Elsevier. 

AMY [152] • The left panels show Raa as a function of Pt for two different centrality bins. They serve as 
proof that both the P^-dependence and the centrality dependence of Raa can be described by all three 
models. Every model has one free parameter that has been fitted: the strong coupling as for AMY, (q) 
or derived parameters for the HT and ASW formalisms. In this particular case a normalization factor 
K = q/{2e^'^) as explained in Eq. ([80]) was fitted for ASW, where e is the local energy density. 

This study confirms the surprisingly large q found in the ASW model compared to other approaches. 
For the case that the quenching strength scales with e^/^ the initial values found for a quark at the 
center of the fireball in central collision are jl52] 

q = 18.5 GeVVfm for ASW , g = 4.5 GeVVfm for HT (109) 

and for the case that the quenching strength scales like the temperature T it is found that 

g = 10 GeVVfm for ASW , g = 2.3 GeVVfm for HT , g = 4.1 GeVVfm for AMY . (110) 

Recall that the rates in AMY are calculated self-consistently as functions of the local temperature so 
there is only one choice to model the space and time dependence. 

This comparison is unique and very valuable since the same initial hard cross sections and the 
same maps for the fireball, from (3+l)-dimensional ideal hydrodynamics were used. Any differences in 
the extracted values of g must be due to differences between the calculations themselves, not due to 
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Figure 31: Comparison between true elliptic flow (solid line) and a 
suggested method to compute reconstructed elliptic flow from data 
vf"^^ (dash-dotted line). Squares represent PHOBOS data for mid- 
central collisions [2841 EHHl [286] . See Ref. [189] for detail. Figure 
reprinted from |189] with permission from the American Physical 
Society. 



differences in implementation. One of the conclusions is that our current knowledge applied to Raa 
leaves a rather large uncertainty in the determination of q. 

The right panel of Fig. [2^ shows Raa as a function of the angle with respect to the reaction plane 
normalized by the average Raa- Due to the large values of q the ASW formalism is more strongly 
dominated by surface emission than the other models. This also leads to a stronger angular modulation 
for non-spherical fireballs. 

4.2 Azimuthal Anisotropics and Elliptic Flow 

Anisotropics in the azimuthal angle 0, especially elliptic flow, contain detailed information on the hot 
and dense QCD bulk matter created at RHIC. We will use the notion of elliptic flow for the second 
harmonic 

[d^cos2^dN/d^dPr 
^'^^"^^ - ^d<pdN/d<pdPT ^ ^ 

for any value of Pt, even if the asymmetry is not produced through hydrodynamic flow. Large elliptic 
flow of the bulk matter produced at RHIC had been observed early on and has led to the claim of 
perfect fluidity of the quark gluon plasma just above Tc [8]. With the advent of viscous relativistic 
hydro codes the interest has shifted toward the goal of quantifying the dissipative transport coefficients, 
in particular the shear viscosity rj. 

In the hydrodynamic regime the asymmetry in the pressure gradient in and out of the reaction plane 
for finite impact parameters drives a larger acceleration in-plane than out-of-plane. Figure [25] shows the 
elliptic fiow V2 as a function of Pt for vr, K and p from an early calculation using ideal hydrodynamics 
by Huovinen et al. |165] . The calculations show remarkable agreement with experimental data. We 
can observe a clear mass ordering: the f 2 for pions is larger than that of kaons which in turn is larger 
that that of protons. The effect is most pronounced at low Pt- This is a natural phenomenon in 
hydrodynamics and comes from the interplay of radial flow in and out of the reaction plane. The 
effects of flow are more pronounced for more massive particles due to their smaller thermal velocities 
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Figure 32: Left panel (a): Centrality dependences of V2 for charged particles at midrapidity in Au+Au 
collisions at ^/SNN = 200 GeV with CGC initial conditions. Open circles (squares) are results with 
(without) eccentricity fluctuations. Right panel (b): The same for Cu+Cu collisions. Figures reprinted 
from |287] with permission from the American Physical Society. 



at a given temperature, which translates into smaller asymmetries. As for spectra the good agreement 
between hydrodynamic models and experimental data is restricted to low Pt as effects of insufficient 
thermalization kick in at higher momenta. Despite the pioneering character of this calculation it is at 
present no longer regarded as being very realistic due to the following issues: (i) the lack of proper 
treatment of chemical vs thermal freeze-out, (ii) the fact that it assumes boost-invariance and ignores 
the dynamics in longitudinal direction, (iii) the absence of final state interactions, (iv) the absence 
of dissipative corrections. As we had pointed out earlier hydrodynamic models that do not take into 
account the difference of chemical and kinetic freeze-outs can not explain the absolute values of proton 
spectra correctly. However this calculation gives us the guidelines for the necessary improvements that 
have been implemented since then. 

Figure [26] shows elliptic flow at RHIC as a function of pseudorapidity r/. This calculation is carried 
out with a (3+l)-dimensional ideal hydro model with two freeze-out processes (labeled PCE = partial 
chemical equilibrium) |202] . These features improve the P^-spectra of protons and give reasonable 
results for f 2(-Pr) for ^r, K and p at midrapidity. However the authors of this study flnd that they can not 
explain the data from RHIC away from midrapidity. Besides generally overestimating transverse flow, 
there are humps in forward and backward rapidity which do not feature in the experiment data |28HI284] . 
It should be noted that the rapidity spectra are described well after two adjustable parameters in the 
initial conditions in longitudinal direction had been flxed. This result might indicate that thermalization 
is reached only very close to midrapidity. 

As in the case of spectra one can also infer from comparison of elliptic flow with data that hydrody- 
namic models work only at < 2 GeV/ c, see Fig. |23 Naturally we expect a transition to a region that 
is dominated by perturbative production, with a recombination region in between. Note that the mech- 
anism of f 2 generation is different for hard probes. An azimuthal asymmetry develops simply because 
of the different amount of material jets have to go through in- and out-of-plane, with smaller opacity 
in-plane. Fig. [271 shows V2 as a function of Pt for vr", K~ , p and charged hadrons for the hydro -|- jet 
model in |267j . In experimental data we observe that the shape of the P^-dependence of V2 changes 
from that of a pure hydro model by bending over at larger Pp, saturating for a while at intermediate 
Pt, and then dipping down at even larger values of Pt- The transition points, as for spectra, depend 
on particle species. In addition, the saturation levels at intermediate Pt shows a peculiar universality 
with baryons and mesons lining up at two different values of V2 which scale like 3:2 [U [5]. 

Hydro + jet models get some of these basic features, but usually can not explain the systematics of 
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Figure 33: Left panel: Integrated V2 as a function of A^part for charged particles in Au+Au collisions 
at y/SNN = 200 GeV, compared to viscous hydrodynamics with shear viscosity for various viscosity- 
to-entropy ratios rj/s and compared to data from PHOBOS |289j . Right panel: V2 as a function of Pt 
for the same system compared to data from STAR [290J. Figures reprinted from [28BJ with permission 
from the American Physical Society. 



transition points and the universal baryons vs meson saturation levels. This is a strong argument for 
the presence of quark recombination at intermediate Pt- Figure [28] shows f 2 as a function of Pt for 
identified particles from the recombination plus fragmentation model advocated in Ref. |251j together 
with early experimental data. Recombination models naturally deliver the universal behavior within 
the baryon and meson groups and indeed predict the valence quark number scaling of f 2 as discussed 
earlier in this review. The data clearly support this stunning feature and have supported the claim that 
traces of collectivity at the parton level can be seen |251j . At higher Pt perturbative hadron production 
takes over. It is generally expected that V2 in the perturbative domain does not depend too much on 
hadron species, but reliable data above 6 GeV/c is scarce. Large differences between hadrons could 
be a sign of hadron- instead of parton-based jet quenching, or they could indicate rapid changes in 
jet chemistry inside a quark gluon plasma |139] . We discuss more theoretical results for V2 at high Pt 
below. Figure [28] shows general agreement of the calculations with data except for the small P^-region 
where the mass splitting is not resolved since no genuine hydrodynamics phase with correct hadron 
masses was used. 

While recombination models lend a helpful hand to hydrodynamics to extend the bulk properties to 
larger Pt, hadronic transport is able to fix our failing understanding of the pseudorapidity dependence 
discussed earlier, and of the centrality scaling of elliptic fiow |191] I211j . In Fig. [29] we show the Pt- 
integrated elliptic fiow for charged hadrons at midrapidity as a function of the number of participating 
nucleons Ai'part- The plot compares results obtained with an ideal (3-|-l)-dimensional hydrodynamic 
model with and without the hadronic cascade model JAM attached as an afterburner from Hirano et 
al. |211] . In this study the authors also compare two different initial conditions for their hybrid model: 
a Glauber model (BGK) and Color Glass Condensate-based model (CGC). We refer the reader to Ref. 
[211] for more details. We observe that the hybrid hydro+JAM model with Glauber initial conditions 
gives the best description of the centrality dependence of elliptic flow. The hydro+JAM model with 
CGC initial conditions which is generally considered more realistic is consistent with experimental data 
only at large A^part- 

The pseudorapidity dependence of V2 for charged hadrons in central, semi-central and peripheral 
Au+Au collisions at ^/SNN = 200 GeV is shown in Fig. [30] By comparing results from pure ideal 
hydrodynamics at decoupling temperatures of T^ec = 100 MeV and T^ec = 169 MeV respectively, we see 
that the bumps at forward and backward rapidities that were already observed in the previous study 
|202] are larger if the hadronic matter is allowed to evolve without dissipation. On the other hand, 
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Figure 34: Estimated band of rj/s form current viscous hydro- 
dynamics models in Tab. |3] and from several theoretical estimates 
discussed in the text. 



V2{f]) from the hybrid hydro plus hadronic cascade model does not allow these structures to build up 
in V2{r]). We can conclude from all of these hybrid studies with hadronic cascades that effects of shear 
viscosity and dissipation in the hadronic phase, and proper final state interactions are not negligible. 
We also note that once more the initial conditions based on the Glauber model have the upper hand over 
CGC initial conditions in comparison with data. Color glass initial conditions would require additional 
dissipation during the QGP phase. However, even this much improved investigation did not take into 
account the effect of event-by-event fluctuations of the geometric shape of the density of the initial 
condition which affects the elliptic flow. 

The importance of including event-by-event fluctuations was discussed systematically by Andrade et 
al. using the hydrodynamic code NEXSPHERIO [189j which includes a Monte-Carlo generator for initial 
conditions. They show that the assumption of symmetry of the particle distribution in relation to the 
reaction plane leads to disagreement between the true and reconstructed elliptic flows and emphasize 
that it is important to have a precise experimental determination of elliptic flow. Their calculated V2 
as a function of pseudorapidity shows very nice agreement with experimental data as shown in Fig. [3T1 
However, they did not connect their hydrodynamic model to a transport code for the hadronic phase. 

The effect of eccentricity fluctuations on the elliptic flow at midrapidity in Au +Au and Cu + 
Cu collisions was recently also investigated in Ref. [287]. Those authors include the effect of initial 
eccentricity fluctuations originating from the nucleon position inside the colliding nuclei both for the 
Glauber model and CGC initial conditions. The effect of eccentricity fluctuations is not very large in 
semi-central Au + Au collisions and it does not shift the values of f 2 closer to experimental data in that 
region. On the other hand, it enhances V2 in Cu + Cu collisions where fluctuations are more important 
because of the smaller system size. As a result V2{i]) from CGC initial conditions with fluctuations can 
describe the experimental data quite well. 

Finally we show results for elliptic flow from one of the recently developed (2+l)-dimensional rela- 
tivistic viscous hydrodynamic codes. Since both the formulation and the implementation of relativistic, 
second order viscous hydrodynamics is non-trivial, qualitative comparison with experimental data is 
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Figure 35: Left panel: V2 as a function of A^part- Right panel: V2 as a function of Pt 
for two different centrality bins. As in Fig. [22] dashed lines refer to reweighted and solid 
lines to non-reweighted ASW quenching weights. STAR and PHENIX data from references 
[282[ I293[ I294j . Figures reprinted from p. 23] with kind permission from Springer Science + 
Business Media. 



just starting. Figure [33] shows results for t>2 of charged particles in Au+Au collisions obtained by Ro- 
matschke and Romatschke |288j together with experimental data from STAR and PHOBOS. This study 
finds that both integrated V2 as a function of centrality, and differential f 2 as a function of Pt need very 
small values of rj/s, 0.03 . . . 0.08, partially violating the conjectured KSS bound. It is too early to reach 
conclusive results in light of the many small improvements that ideal hydrodynamics needed to imple- 
ment over a decade to become a reliable tool. E.g. the study mentioned here considers only one type of 
initial condition (Glauber) and it does not implement final state interactions |288] FI Nevertheless, the 
first efforts have jump started the era of qualitative studies using dissipative hydrodynamics |178l I291j . 

In Fig. [34] we show a band for favored values of the shear viscosity over entropy ratio rj/s which is 
extracted from the recent viscous hydrodynamic models listed in Tab. [3] We also show results from 
lattice QCD [225], UrQMD [231], pQCD (Fig. 10 in Ref. [178]) and phenomenological analyses [292] for 
reference. This rj/s band from viscous hydrodynamics should not be considered as a conclusive result, 
but it highlights an interesting preliminary finding: The rj/s extracted from comparison of viscous 
hydrodynamics with RHIC data is generally in the vicinity of the KSS bound. From the comparison 
with UrQMD and lattice QCE0, we would conclude that those small values of rj/s must come from the 
QGP phase. These facts support the hypothesis of a sQGP at RHIC. In these viscous hydrodynamic 
calculations the temperature dependence of viscosities is usually not taken into account, as shown in 
Fig. [3U The temperature dependence of transport coefficients may be not very significant in the QGP 
phase at RHIC. However the shear viscosity of the hadron phase seems to be much larger than that of 
the QGP phase which suggests that it is necessary to take this difference into account when discussing 
phenomena related to the QCD phase transition |167j . 

Let us come back to the problem of azimuthal anisotropy at large Pt- We have already discussed the 
basic mechanism how perturbative hadron production with final state interaction in a medium can lead 
to positive t>2. We want to close this subsection by showing the V2 obtained in the study by Dainese, 

^ Luzum and Romatschke investigate the initial and freeze-out temperature dependence of eUiptic flow, using a 
Glauber-based model and CGC, see Sec. 13.2.21 

''The calculation is performed with SU(3) pure gauge theory. 
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Figure 36: HBT radii for negative pions in (from 
the top) side, out, and long directions and the ratio 
-RoutZ-Rside (bottom) as function of pair momentum 
Kt- Sohd and dashed hues correspond to hydro cal- 
culations using partial chemical equilibrium (PCE) 
and full chemical equilibrium (CE) resp. STAR and 
PHENIX data from Refs. [2961 1297] Figure reprinted 
from |202j with permission from the American Phys- 
ical Society. 
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Figure 37: HBT radii and the ratio 
-RoutZ-Rside for four different scenarios for ini- 
tial condition and freeze-out, see text for 
details. Data for pions from STAR and 
PHENIX are shown [2961 IM]- Figure 
reprinted from |215j with permission from 
the American Physical Society. 



Loizides and Paic, using the ASW formalism, that we already discussed for their results on Raa in the 
previous subsection |123j . Fig. [551 shows V2 as a function of centrahty (left panel) and as a function of 
Pt (right panel) with data from STAR and PHENIX. We notice that the asymmetry from the difference 
in opacity in- and out-of-plane can be sizable, but the calculation still underestimates most of the data. 
The large V2 measured by experiments at large Pt has long been puzzling, but the experimental data 
also exhibits large error bars. Note that q was fixed in order to describe Raa which leaves no free 
parameter in this study. 

We want to remind the reader of the right panel in Fig. [2H In that study Raa was investigated as 
a function of the azimuthal angle (f) with respect to the reaction plane, not just integrated over 0. In 
principle Raa{4>, Pt) has more differential information than either V2 or integrated Raa- We recall that 
the ASW formalism exhibited the strongest angular modulation and has therefore the largest f 2 in that 
comparative study. 
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Figure 38: Gaussian HBT radii -Rout, -Rside, and -Riong from several hydrody- 
namic calculations explained in the text, together with data from STAR (red 
stars). Figure reprinted from [298] with permission from the American Physi- 
cal Society. 

4.3 Two-Particle Correlations 

Before RHIC started up two-pion Hanbury Brown-Twiss (HBT) interferometry was believed to give us 
a clear signature of the QCD phase transition. We expected that an enhancement of the ratio of the 
inverse width of the pion correlation function in out-direction to that in side-direction, which results 
from a prolonged life time of the fire ball with a phase transition, could be observed at RHIC [295j . 
However what we find in experimental data from RHIC are HBT radii that are almost the same as 
those measured at SPS [280] . Furthermore, most present hydrodynamic models can not describe the 
measured HBT radii correctly although they fit both spectra and elliptic flow. 

Figure [36] shows the HBT radii -Rside, -Rout, -Riong and the ratio -RoutZ-Rside for negative pions calcu- 
lated from (3+l)-dimensional relativistic hydrodynamics for both partial chemical equilibrium (PCE) 
and chemical equilibrium (CE) in the hadronic phase by Hirano et al. |202j . Except for -Riong the 
hydrodynamic calculations fail to reproduce the experimental data quantitatively. -Rside (-Rout) from 
the hydrodynamic model underestimates (overestimates) the experimental data. This leads to a large 
discrepancy in the ratio of -Rout /-Rside between hydrodynamic models and experimental data, which is 
called the HBT puzzle. The same tendency can be seen in other hydrodynamic calculations. At face 
value this means that the expansion of the fireball happens more rapidly in a shorter time than in 
hydrodynamic models. Partial chemical equilibrium pushes the hydrodynamic calculations closer to 
experimental data, however it is not enough to solve the HBT puzzle. 

Two improvements are taken into account in Ref . [215] . For one, event-by-event fluctuations in the 
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Figure 39: Iaa for charged hadrons as a function of A'^part with trigger particle having 4 < 
Pt < 6 GeV/c compared to data from STAR |299j. Figure reprinted from |123j with kind 
permission from Springer Science and Business Media. 



initial conditions, and secondly continuous emission instead of the sudden freeze-out process which is 
usually used in hydrodynamic models. The results are shown in Fig. [371 The same radii and ratios 
as before are shown in four cases: sudden freeze-out with averaged initial condition (FOl), sudden 
freeze-out with fluctuating initial condition (F02), continuous emission with averaged initial conditions 
(CEl) and continuous emission with fluctuating initial conditions (CE2). The realistic treatments of 
initial conditions and freeze-out bring a significant improvement in -Rout, which also turns into a better 
agreement of Rout/RsiAe- However, small discrepancies between hydrodynamic calculations and the 
experimental data on Rout remain. 

A solution to the HBT puzzle was finally proposed by Pratt |298j . He suggests that the discrepancies 
are not from a single shortcoming of hydrodynamic calculations, but from a combination of several 
effects: mainly prethermalized acceleration, equations of state with inadequate stiffness, and lacking 
viscosity. Figure [38] shows results from his calculations for the HBT radii together with data from 
STAR (red stars). The results of many calculations are shown with the extremes being from a hydro 
calculation with a first-order phase transition without pre-thermal flow and without viscosity (black 
squares and line), and the gradual improvements culminate in a calculation with stiffer equation of 
state and pre-thermal flow and viscosities included (black circles and line). 

At high Pt 2-particle correlations are measured in heavy ion collisions not for their information 
on quantum interference but in order to establish kinematic correlations. Di-hadron measurements at 
RHIC have made it possible to study some properties of jets even though true jet reconstruction is 
remarkably difficult. Correlation measurements can be recorded by using pairs T — A where T is a 
trigger hadron and A the associated particle. Most studies at intermediate and high Pt have been 
carried out using such triggered correlations, looking at the per-trigger yield of associated particles 

as a function of the relative azimuthal angle A0 between trigger and associate. Here Pt and Pa are, in 
deviation from our usual notation, the transverse momenta of the trigger and associated hadrons, resp. 
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Figure 40: x^"fits of the parameter K in the ASW model of energy loss from Raa and Iaa 
for three different extrapolations of jet quenching into the pre-equilibrium phase. Figures 
reprinted from |151] with permission from lOP. 



One can then proceed and define a nuclear modification factor for the associated yield 



i-AA 



Ypp{Pt,Pa,A(P) 



(113) 



Like Raa we expect Iaa to be unity for a loose superposition of proton-proton collisions. Another 
potentially useful observable are triggered fragmentation functions. They can be derived from per- 
trigger yields by rewriting the dependence as one on a "momentum fraction" z = Pa/Pt and integrating 
over a narrow bin around A(f) = n. 

The amount of data collected on this topic would merit its own review article. We will focus on a 
few selected examples, mainly to achieve our goal to better constrain the transport coefficient q. We 
start by showing the result from Ref. [123] which is now well-known from our discussions on Raa and 
V2- Fig. [39] shows their result for Iaa as a function of centrality in the ASW model of energy loss 
compared to data from STAR. The data fall well into the large error band given by the uncertainty 
from the reweighted vs the non-reweighted quenching weights. 

Let us recall the right panel of Fig. |23l In that study the Higher Twist formalism was used together 
with a calculation of hard processes at NLO accuracy to extract the stopping power eo. As shown in 
the figure, fitting both Iaa and Raa leads to consistent values for eo. This is somewhat different for the 
recent study by Armesto et al. in the ASW model where the consistency of q extracted from Raa and 
Iaa data is very sensitive to details of the modeling |151] . Fig. HO] shows the result for K = g/(2e^/^) 
extracted from Raa and Iaa with varying treatments of the quenching during the pre-equilibrium phase. 
The three cases shown are (i) no quenching before the QGP formation time Tq, (ii) constant q for t < tq 
and (iii) q increasing as t~^^^ toward r = 0. The results is a reminder that the determination of q from 
data has larger uncertainties. 



4.4 Photons 

Electromagnetic probes are important means to obtain information from the quark gluon plasma. Pho- 
tons and lepton pairs carry the information of the whole time-evolution of heavy ion collisions: initial 
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Figure 41: Direct photon production in 
Au+Au collisions at a/s/vjv = 200 GeV 
for five different centrality bins compared 
to data from PHENIX [M]. Figure 
reprinted from |300j with permission from 
the American Physical Society. 
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Figure 42: Nuclear modification factor Raa of direct 
photons as in Fig. SH Left panel: three different cen- 
trality bins are shown. Right panel: calculations with 
and without energy loss of partons. Figures reprinted 
from [300j with permission from the American Physical 
Society. 



hard collisions, thermalization, expansion, hadronization, and freeze-out. We focus here on photons 
which encode detailed information of every process which occurs in heavy ion collisions. A schematic 
account of photon sources has already been given, but we want to list the sources of direct photons 
once more in one place: (i) prompt hard photons from initial collisions, (ii) vacuum bremsstrahlung 
(fragmentation) photons, (iii) photons from jet conversions and medium-induced bremsstrahlung (iv) 
thermal photons from quark gluon plasma (and the hot hadronic phase). This list also orders photons 
according to which momentum regime they are most important for (from the largest Pt to the smallest). 

We can deduce the temperature of the fireball from thermal radiation (and check our understanding 
of the fireball evolution), and we can infer information about the medium density and q from jet 
conversions and induced photon bremsstrahlung. Photons have also great importance as triggers in 
photon-hadron correlation studies at large Pt- Triggered fragmentation functions with photon triggers 
come close to real fragmentation functions since the photon, if the source is dominated by initial hard 
photons, carries the same transverse momentum as its initial partner parton jTl4] . This opens a way to 
measure real medium-modified fragmentation functions, however one has to be cautious because of the 
many other sources that weaken this kinematic link between the photon and the jet on the other side. 

Figures SU and Il2] show the centrality dependence and the nuclear modification factor of direct 
photons in Au+Au collisions at ^Js^n = 200 GeV, respectively, in the study by Liu et al. [300] . They 
include all four sources mentioned above and compute thermal radiation by using (3+l)-dimensional 
hydrodynamic fireball evolution. Their calculation describes the direct photon data from PHENIX very 
well. The same can be said about the calculation of the McGill group |302] . They use all four sources 
of photons as well, with energy loss of jets and induced photon bremsstrahlung taken care of by the 
AMY formalism. A successful fit of hadron and photon data together is hence a crucial consistency test 
for the AMY formalism. Fig. |13] shows their result for both single inclusive photon spectra and Raa 
for large and intermediate Pt (thermal radiation is therefore omitted in this case). 
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Figure 43: Left panel: direct photon spectra for central Au+Au collisions at RHIC using 
contributions from hard collisions, bremsstrahlung and jet conversions together with data 
from PHENIX [303] . Right panel: Raa for direct photons from the same sources. Figures 
reprinted from j3U2] with permission from the American Physical Society. 



Fig. HH shows a calculation of the McGill group including thermal photons according to Arnold, 
Moore, and Yaffe [304] . In that study the authors also calculated the azimuthal asymmetry V2 of 
direct photons. The V2 of thermal photons and of fragmentation photons is expected to be positive, 
reminiscent of the V2 of bulk hadrons and high-P^ hadrons respectively. However, photons from jet 
conversions and induced bremsstrahlung add a negative contribution to the direct photon f 2 as they are 
produced more abundantly in the direction where the medium is thicker [305]. The left panel of Fig. 
HH shows different contributions to V2 with the contribution from jet conversions indeed being negative. 
Experimental data from PHENIX has been inconclusive so far due to large error bars. Measurements 
of negative values of V2 would be direct evidence for jet conversions. 

Fig. |45] shows photon-triggered fragmentation functions Daa obtained from photon-hadron correla- 
tions calculated in Ref. [302] with the AMY formalism together with data from PHENIX. This is one of 
the emerging examples of jet tomography with photon triggers as originally envisioned in |114] . In this 
study the data is described reasonably well with the parameter in the AMY energy loss (the coupling 
constant ag) set to one consistent value for all observables. 

5 Summary and Conclusions 

Let us briefly summarize. We have laid out the foundations of perturbative QCD and how they can be 
used to understand data from the Relativistic Heavy Ion Collider. We explained the concept of collinear 
factorization widely used in collider physics, and modifications when applied to colliding nuclei. We 
have then focused on final state interactions, and in particular on leading particle energy loss. Along 
the way have discussed several popular model calculations for energy loss. We have also pointed out 
some of their basic assumptions in which they differ. All of them have apparent shortcomings, and 
even if jet quenching at RHIC is perfectly perturbative, we should not expect all of these calculations 
to apply. However, all of them explain jet quenching well on a qualitative level. This can be seen 
positive, we are on the right track, or negative, we are not in a situation where we can confidently 
falsify one or all of those calculations, due to a long list of uncertainties. We are also not quite ready 
yet to exclude non-perturbative quenching scenarios from data. In fact, it is very well possible that a 
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Figure 44: Left panel: Direct photon spectra from all known sources for central Au+Au 
collisions using the AMY formalism consistently. Right panel: direct photon V2 as function 
of Pt- Figures reprinted from |3U4] with permission from the American Physical Society. 



mixture of perturbative and non-perturbative quenching co-exist in different temperature ranges. 

We showed that all models fit the data on single hadron suppression after adjusting a single param- 
eter, including the systematics of the Pt- and centrality dependence of Raa- On the other hand the 
quenching strengths extracted from the different models is quite different. The maximum values of q in 
central collisions that have been found lie in the range ~ 2 ... 15 GeV^/fm. Other observables like V2 
at high Pt or I^a have note yet cut down this range as it has become clear that the results are quite 
sensitive, e.g. to details of the treatment of the underlying fireball. 

There is a focused effort to attack the open questions within the framework of the TECHQM 
collaboration and the newly founded JET collaboration. In the future we need rigorous comparisons 
and a vetting process of different calculations, a realistic modeling of the fireball and realistic assessments 
of uncertainties in order to start to falsify certain assumptions and to reduce the error bar on q. We also 
have to include more realistic options, like combining perturbative quenching with final state effects in 
the hadronic phase. 

We have also reviewed the successful story of hydrodynamics at RHIC, and the novel developments 
of the past few years that have brought the first relativistic viscous hydro codes and hybrid models. 
We have discussed a long list of observables that require the presence of a thermalized and collectively 
moving quark gluon plasma phase, that can be beautifully explained by hydrodynamic calculations and 
their extensions. The most convincing single observations are the mass ordering of flow and elliptic 
flow, and the large size of elliptic flow at RHIC, which are suggested by hydrodynamics and observed 
in the data. Starting from basic concepts of ideal and viscous hydrodynamics we have moved along a 
variety of steps that brought vast improvements to our understanding of hydrodynamics, including the 
arrival of 3+1 dimensional codes, the correct treatment of freeze-out, etc. 

Like in the case of pQCD and jet quenching, hydrodynamics has had some difficulties entering 
a period of precision measurements. Currently, the equation of state from the comparison between 
hydrodynamic analyses and experimental data favor a first order deconfinement transition, while lattice 
QCD prefers a crossover transition. Clearly, more systematic phenomenological studies, are needed 
that try to start from more relaxed assumptions, e.g. by including initial radial fiow and dissipative 
corrections. It is not until both sides, phenomenological analyses of data and theoretical tools like 
lattice QCD, show the same conclusion that we can pin down the QCD phase diagram. Routine runs 
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Figure 45: Photon-triggered fragmentation functions Dj^a for three different trigger windows 
in the same model as in Fig. US] compare to data from PHENIX [306l| . Figure reprinted from 
[302] with permission from the American Physical Society. 



of fully 3+1- dimensional, viscous hydrodynamics are around the corner and they should bring us a step 
closer towards an understanding of QCD thermodynamics and transport coefficients. Current estimates 
of the shear viscosity range from the 1 to about 4 times the KSS bound of s/Atc. 

One of the big challenges in this field is the huge amount of data from RHIC that has remained 
almost untouched by comprehensive theoretical calculations. There are many results on dihadron 
azimuthal correlations at high Pt |307j . three-particle azimuthal correlations, the ridge structure |308j . 
dihadrons with respect to the reaction plane |309j . and the emerging full-jet reconstruction |310l 1311] . 
All of these results should fit into a tent anchored on both sides by pQCD and hydrodynamics. Multi- 
module modeling, carefully tested and endowed with realistic estimates of theoretical uncertainties 
should eventually be able to explain those features and in the process the data should narrow down the 
error bars on transport coefficients significantly. 

On the experimental side two new challenges will arrive shortly. Energy scans will be conducted at 
RHIC |312] and supplemented by SPS results and the program at the new FAIR and NICA facilities 
[313] . These programs will provide new data on the QCD phase diagram at high net baryon densities. 
Reliable hydro+cascade hybrid models will be given a chance to look for the QCD critical point and 
a change in the nature of the phase transition. On the other hand the LHC will vastly improve the 
reach of hard probes by colliding ions at unprecedented center of mass energies. Reconstructed jets and 
abundant high-P^ probes will provide much stricter tests of our understanding of how jets interact with 
quark gluon plasma. 

As the first 10 years of the high energy heavy ion era come to an end we find that we have gained 
much knowledge. It prepares us for the next decade in which the QCD phase diagram and the quark 
gluon plasma phase will be mapped out quantitatively. 
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